Input

Read in individual subject rating data.
Columns are:
id (subject ID, possibly recode for public release)
Age (in years)
Sex (Male or Female)
InjectionUser (TRUE for participants who report a history of injectin methamphetamine, opioids, or both)

ImageSet: control, opioid, or meth
File: file name
Category: subcategory within ImageSet. Possible categories are:

ImageSetFile:

Order: the order the images were seen for this subject, from 0 to 359, spans across two visits.
Visit: V1 or V2 for the two visits
time: the order the images were seen for an individual visit, from 0 to 179
RatingType (which question this response is for: valence, arousal, craving, related, typicality)
Questions were worded:

Value: the value for this rating

#rating_data <- read.csv('OpioidMethRatingSubjectData.csv')
#clinical_data <- read.csv('DCR_ClinicalData.csv')
rating_data <- read.csv('DCR_RatingSubjectData-anonymized.csv')
clinical_data <- read.csv('DCR_Clinical-anonymized.csv')
subject_data <- merge(rating_data, clinical_data)

Individual Subject Plots

library(ggplot2)
library(gridExtra)
library(grid)
library(ggpubr)
## Warning: package 'ggpubr' was built under R version 3.5.3
## Loading required package: magrittr
#put levels in the order we like
subject_data$RelatedCategory <- factor(subject_data$RelatedCategory, levels = c('Neither', 'Meth', 'Opioid', 'Both'))

#subject_data$Category[subject_data$Category == 'meth_injection_instrument'] <- 'meth_instrument'
#cols <- c('control' = 'green', 'meth' = 'red', 'opioid' = 'blue', 'Control' = 'green', 'Meth' = 'red', 'Opioid' = 'blue')
cols <- c('control' = '#7CAE00', 'meth' = '#F8766D', 'opioid' = '#00BFC4', 'Control' = '#7CAE00', 'Meth' = '#F8766D', 'Opioid' = '#00BFC4')

for (subject in unique(subject_data$id)){
  data_to_plot <- subject_data[subject_data$id == subject,]
  
  
  p1 <- ggplot(data = data_to_plot[data_to_plot$RatingType == 'valence',]) + 
    geom_point(aes(x = Order, y = Value, color = ImageSet)) + ggtitle('Valence') +
    xlab('') + ylab('') + theme(plot.title = element_text(hjust = 0.5)) +
    scale_x_continuous(breaks = c(0, 90, 180, 270, 360), limits = c(0, 360)) + 
    scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9)) +
    theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
    scale_color_manual(values = cols)
    
  p2 <- ggplot(data = data_to_plot[data_to_plot$RatingType == 'arousal',]) + 
    geom_point(aes(x = Order, y = Value, color = ImageSet)) + ggtitle('Arousal') +
    xlab('') + ylab('') + theme(plot.title = element_text(hjust = 0.5)) +
    scale_x_continuous(breaks = c(0, 90, 180, 270, 360), limits = c(0, 360)) + 
    scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9)) +
    theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
    scale_color_manual(values = cols)
  
  p3 <- ggplot(data = data_to_plot[data_to_plot$RatingType == 'craving',]) + 
    geom_point(aes(x = Order, y = Value, color = ImageSet)) + ggtitle('Craving') +
    xlab('') + ylab('') + theme(plot.title = element_text(hjust = 0.5)) +
    scale_x_continuous(breaks = c(0, 90, 180, 270, 360), limits = c(0, 360)) + 
    scale_y_continuous(breaks = c(0, 25, 50, 75, 100), limits = c(0, 100)) +
    theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
    scale_color_manual(values = cols)
  
  p4 <- ggplot(data = data_to_plot[data_to_plot$RatingType == 'typicality',]) + 
    geom_point(aes(x = Order, y = Value, color = ImageSet)) + ggtitle('Typicality') +
    xlab('') + ylab('') + theme(plot.title = element_text(hjust = 0.5)) +
    scale_x_continuous(breaks = c(0, 90, 180, 270, 360), limits = c(0, 360)) + 
    scale_y_continuous(breaks = c(0, 25, 50, 75, 100), limits = c(0, 100)) +
    theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
    scale_color_manual(values = cols)
  
  p5 <- ggplot(data = data_to_plot[data_to_plot$RatingType == 'related',]) + 
    geom_point(aes(x = Order, y = RelatedCategory, color = ImageSet)) + ggtitle('Relatedness') +
    xlab('') + ylab('') + theme(plot.title = element_text(hjust = 0.5)) +
    scale_x_continuous(breaks = c(0, 90, 180, 270, 360), limits = c(0, 360)) + 
    theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
    scale_color_manual(values = cols)
  
  filename <- paste(subject, 'responses.png', sep = '_')
  png(filename, width = 1400, height = 400)
  print(annotate_figure(ggarrange(p5, p1, p2, p3, p4, ncol = 5, nrow = 1, common.legend = TRUE, legend = 'right'), top = text_grob(subject, color = 'red')))
  dev.off()
  print(p1 + ggtitle(subject))
  print(p2 + ggtitle(subject))
  print(p3 + ggtitle(subject))
  print(p4 + ggtitle(subject))
  print(p5 + ggtitle(subject))
   
}

Summary Function

library(reshape2)

one_summary <- function(dataset, prefix, grouping_field){
  #dataset will be filtered down to either: all subjects, injection users, non-injection users
  #prefix will be added to the output variable names, for example: allsubjects_arousal_mean or injectionusers_arousal_mean
  
  #get raw mean for each image
  mean_values <- aggregate(dataset[, c('Value')], 
                            by = list(RatingType = dataset$RatingType, Group = dataset[, grouping_field]), 
                            FUN = mean, na.rm = TRUE)
  wide_mean_values <- dcast(mean_values, formula = Group ~ RatingType, value.var = 'x')

  
  category_data <- dataset[dataset$RatingType == 'related',]
  #so we can easily get a count of 'Neither', 'Both', 'Opioid', 'Meth'
  category_data$count_col <- 1
  #get counts of 'neither' ratings
  counts_neither <- aggregate(category_data[category_data$RelatedCategory == 'Neither', c('count_col')], 
                            by = list(Group = category_data[category_data$RelatedCategory == 'Neither', grouping_field]), 
                            FUN = sum, na.rm = TRUE)
  names(counts_neither) <- c('Group', 'CountNeither')
  #'get count of 'meth' ratings
  counts_meth <- aggregate(category_data[category_data$RelatedCategory == 'Meth', c('count_col')], 
                            by = list(Group = category_data[category_data$RelatedCategory == 'Meth', grouping_field]), 
                            FUN = sum, na.rm = TRUE)
  names(counts_meth) <- c('Group', 'CountMeth')
  
  all_counts <- merge(counts_neither, counts_meth, all = TRUE)
  #get count of 'opioid' ratings
  counts_opioid <- aggregate(category_data[category_data$RelatedCategory == 'Opioid', c('count_col')], 
                            by = list(Group = category_data[category_data$RelatedCategory == 'Opioid', grouping_field]), 
                            FUN = sum, na.rm = TRUE)
  names(counts_opioid) <- c('Group', 'CountOpioid')
  
  all_counts <- merge(all_counts, counts_opioid, all = TRUE)
  #get count of 'both' ratings
  counts_both <- aggregate(category_data[category_data$RelatedCategory == 'Both', c('count_col')], 
                            by = list(Group = category_data[category_data$RelatedCategory == 'Both', grouping_field]), 
                            FUN = sum, na.rm = TRUE)
  names(counts_both) <- c('Group', 'CountBoth')
  
  all_counts <- merge(all_counts, counts_both, all = TRUE)
  #any 'NAs that appear here are ratings that never happened--e.g. control images that were never rated as related to meth
  all_counts[is.na(all_counts)] <- 0
  count_cols <- c('CountNeither', 'CountMeth', 'CountOpioid', 'CountBoth')
  all_counts$FractionNeither <- all_counts$CountNeither / rowSums(all_counts[, count_cols])
  all_counts$FractionMeth <- all_counts$CountMeth / rowSums(all_counts[, count_cols])
  all_counts$FractionOpioid <- all_counts$CountOpioid / rowSums(all_counts[, count_cols])
  all_counts$FractionBoth <- all_counts$CountBoth / rowSums(all_counts[, count_cols])
  
  #this will be a number from -1 to 1: -1 for meth, 1 for opioid, and 0 for meth=opioid and many both/neither
  all_counts$MethToOpioid <- (all_counts$CountOpioid - all_counts$CountMeth) / (all_counts$CountMeth + all_counts$CountOpioid + all_counts$CountBoth + all_counts$CountNeither)
  
  
  #add prefix to all_counts column names
  names(all_counts)[!(names(all_counts) %in% c('Group'))] <-
    paste(prefix, names(all_counts)[!(names(all_counts) %in% c('Group'))], sep = '_')
  
  #add _mean to the column names
  names(wide_mean_values)[names(wide_mean_values) %in% c('arousal', 'craving', 'related', 'typicality', 'valence', 'valence_fromneutral')] <- 
    paste(prefix, names(wide_mean_values)[names(wide_mean_values) %in% c('arousal', 'craving', 'related', 'typicality', 'valence', 'valence_fromneutral')], 'mean', sep = '_')
  
  #get raw SD for each image
  sd_values <- aggregate(dataset[, c('Value')], 
                              by = list(RatingType = dataset$RatingType, Group = dataset[, grouping_field]), 
                              FUN = sd, na.rm = TRUE)
  wide_sd_values <- dcast(sd_values, formula = Group ~ RatingType, value.var = 'x')

  #add _sd to the column names
  names(wide_sd_values)[names(wide_sd_values) %in% c('arousal', 'craving', 'related', 'typicality', 'valence', 'valence_fromneutral')] <- 
    paste(prefix, names(wide_sd_values)[names(wide_sd_values) %in% c('arousal', 'craving', 'related', 'typicality', 'valence', 'valence_fromneutral')], 'sd', sep = '_')


  #merge them into one table
  merged_summary_values <- merge(wide_mean_values, wide_sd_values)

  merged_summary_values <- merge(merged_summary_values, all_counts)
  
  return(merged_summary_values)
}

Simple Means/SDs for each image

Will report means/SDs for each image among all subjects, injection users, and non-injection users.

#also get abs(valence - 5) for each rating--distance from neutral
valence_rows <- subject_data[subject_data$RatingType == 'valence',]
valence_rows$RatingType <- 'valence_fromneutral'
valence_rows$Value <- abs(valence_rows$Value - 5)
subject_data <- rbind(subject_data, valence_rows)


all_data_summaries <- one_summary(subject_data, 'allsubjects', 'ImageSetFile')

injection_user_summaries <- one_summary(subject_data[subject_data$InjectionUser,], 'injectionusers', 'ImageSetFile')
noninjection_user_summaries <- one_summary(subject_data[!subject_data$InjectionUser,], 'noninjectionusers', 'ImageSetFile')

merged_summary_values <- merge(all_data_summaries, injection_user_summaries)
merged_summary_values <- merge(merged_summary_values, noninjection_user_summaries)

#now, add in HSV values computed by MATLAB
hsv_values <- read.csv('DCR_hsv_values.csv')

hsv_values$Group <- paste(hsv_values$ImageSet, hsv_values$File)
merged_summary_values <- merge(merged_summary_values, hsv_values)

#add in category column
category_data <- subject_data[!duplicated(subject_data$ImageSetFile), c('ImageSetFile', 'Category')]
names(category_data) <- c('Group', 'Category')

merged_summary_values <- merge(merged_summary_values, category_data)

merged_summary_values$p_cravingdelta_injectionVnon <- NA
###add in p-value for difference in craving between injection and non-injection users for each image
for (image in unique(subject_data$ImageSetFile)){
  injection_craving_ratings <- subject_data[(subject_data$ImageSetFile == image) & subject_data$InjectionUser & subject_data$RatingType == 'craving', 'Value']
  noninjection_craving_ratings <- subject_data[(subject_data$ImageSetFile == image) & (!subject_data$InjectionUser) & subject_data$RatingType == 'craving', 'Value']
  this_test <- t.test(injection_craving_ratings, noninjection_craving_ratings)
  merged_summary_values[merged_summary_values$Group == image, 'p_cravingdelta_injectionVnon'] <- this_test$p.value
}

Category means/sds for the same measures as before

###add in category means/SDs for the same measures

all_data_summaries_categories <- one_summary(subject_data, 'allsubjects', 'Category')
injection_user_summaries_categories <- one_summary(subject_data[subject_data$InjectionUser,], 'injectionusers', 'Category')
noninjection_user_summaries_categories <- one_summary(subject_data[!subject_data$InjectionUser,], 'noninjectionusers', 'Category')

merged_summary_values_categories <- merge(all_data_summaries_categories, injection_user_summaries_categories)
merged_summary_values_categories <- merge(merged_summary_values_categories, noninjection_user_summaries_categories)

###add in mean/sd of hue/saturation/value mean/sd
this_mean_values <- aggregate(merged_summary_values[, c('hue_mean')], 
                            by = list(Category = merged_summary_values$Category), 
                            FUN = mean, na.rm = TRUE)
names(this_mean_values) <- c('Group', 'hue_mean')
this_sd_values <- aggregate(merged_summary_values[, c('hue_mean')], 
                            by = list(Category = merged_summary_values$Category), 
                            FUN = sd, na.rm = TRUE)
names(this_sd_values) <- c('Group', 'hue_sd')
merged_summary_values_categories <- merge(merged_summary_values_categories, this_mean_values)
merged_summary_values_categories <- merge(merged_summary_values_categories, this_sd_values)

##saturation
this_mean_values <- aggregate(merged_summary_values[, c('saturation_mean')], 
                            by = list(Category = merged_summary_values$Category), 
                            FUN = mean, na.rm = TRUE)
names(this_mean_values) <- c('Group', 'saturation_mean')
this_sd_values <- aggregate(merged_summary_values[, c('saturation_mean')], 
                            by = list(Category = merged_summary_values$Category), 
                            FUN = sd, na.rm = TRUE)
names(this_sd_values) <- c('Group', 'saturation_sd')
merged_summary_values_categories <- merge(merged_summary_values_categories, this_mean_values)
merged_summary_values_categories <- merge(merged_summary_values_categories, this_sd_values)

##value
this_mean_values <- aggregate(merged_summary_values[, c('value_mean')], 
                            by = list(Category = merged_summary_values$Category), 
                            FUN = mean, na.rm = TRUE)
names(this_mean_values) <- c('Group', 'value_mean')
this_sd_values <- aggregate(merged_summary_values[, c('value_mean')], 
                            by = list(Category = merged_summary_values$Category), 
                            FUN = sd, na.rm = TRUE)
names(this_sd_values) <- c('Group', 'value_sd')
merged_summary_values_categories <- merge(merged_summary_values_categories, this_mean_values)
merged_summary_values_categories <- merge(merged_summary_values_categories, this_sd_values)

#put ImageSet back in
library(stringr)
merged_summary_values_categories$ImageSet <- str_split_fixed(merged_summary_values_categories$Group, '_', 2)[,1]

###add in p-value for difference in craving between injection and non-injection users for each category
merged_summary_values_categories$p_cravingdelta_injectionVnon <- NA
for (category in unique(subject_data$Category)){
  injection_craving_ratings <- subject_data[(subject_data$Category == category) & subject_data$InjectionUser & subject_data$RatingType == 'craving', 'Value']
  noninjection_craving_ratings <- subject_data[(subject_data$Category == category) & (!subject_data$InjectionUser) & subject_data$RatingType == 'craving', 'Value']
  this_test <- t.test(injection_craving_ratings, noninjection_craving_ratings)
  merged_summary_values_categories[merged_summary_values_categories$Group == category, 'p_cravingdelta_injectionVnon'] <- this_test$p.value
}




library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.3
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:gridExtra':
## 
##     combine
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(data.table)
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## The following objects are masked from 'package:reshape2':
## 
##     dcast, melt
merged_summary_values <- rbindlist(list(merged_summary_values, merged_summary_values_categories), fill = TRUE, use.names = TRUE)
merged_summary_values <- data.frame(merged_summary_values)
all_data_summaries_sets <- one_summary(subject_data, 'allsubjects', 'ImageSet')
injection_user_summaries_sets <- one_summary(subject_data[subject_data$InjectionUser,], 'injectionusers', 'ImageSet')
noninjection_user_summaries_sets <- one_summary(subject_data[!subject_data$InjectionUser,], 'noninjectionusers', 'ImageSet')

merged_summary_values_sets <- merge(all_data_summaries_sets, injection_user_summaries_sets)
merged_summary_values_sets <- merge(merged_summary_values_sets, noninjection_user_summaries_sets)


###HSV code


###add in mean/sd of hue/saturation/value mean/sd
this_mean_values <- aggregate(merged_summary_values[, c('hue_mean')], 
                            by = list(ImageSet = merged_summary_values$ImageSet), 
                            FUN = mean, na.rm = TRUE)
names(this_mean_values) <- c('Group', 'hue_mean')
this_sd_values <- aggregate(merged_summary_values[, c('hue_mean')], 
                            by = list(ImageSet = merged_summary_values$ImageSet), 
                            FUN = sd, na.rm = TRUE)
names(this_sd_values) <- c('Group', 'hue_sd')
merged_summary_values_sets <- merge(merged_summary_values_sets, this_mean_values)
merged_summary_values_sets <- merge(merged_summary_values_sets, this_sd_values)

##saturation
this_mean_values <- aggregate(merged_summary_values[, c('saturation_mean')], 
                            by = list(ImageSet = merged_summary_values$ImageSet), 
                            FUN = mean, na.rm = TRUE)
names(this_mean_values) <- c('Group', 'saturation_mean')
this_sd_values <- aggregate(merged_summary_values[, c('saturation_mean')], 
                            by = list(ImageSet = merged_summary_values$ImageSet), 
                            FUN = sd, na.rm = TRUE)
names(this_sd_values) <- c('Group', 'saturation_sd')
merged_summary_values_sets <- merge(merged_summary_values_sets, this_mean_values)
merged_summary_values_sets <- merge(merged_summary_values_sets, this_sd_values)

##value
this_mean_values <- aggregate(merged_summary_values[, c('value_mean')], 
                            by = list(ImageSet = merged_summary_values$ImageSet), 
                            FUN = mean, na.rm = TRUE)
names(this_mean_values) <- c('Group', 'value_mean')
this_sd_values <- aggregate(merged_summary_values[, c('value_mean')], 
                            by = list(ImageSet = merged_summary_values$ImageSet), 
                            FUN = sd, na.rm = TRUE)
names(this_sd_values) <- c('Group', 'value_sd')
merged_summary_values_sets <- merge(merged_summary_values_sets, this_mean_values)
merged_summary_values_sets <- merge(merged_summary_values_sets, this_sd_values)

#put ImageSet back in
library(stringr)
merged_summary_values_sets$ImageSet <- merged_summary_values_sets$Group

###add in p-value for difference in craving between injection and non-injection users for each category
merged_summary_values_sets$p_cravingdelta_injectionVnon <- NA
for (ImageSet in unique(subject_data$ImageSet)){
  injection_craving_ratings <- subject_data[(subject_data$ImageSet == ImageSet) & subject_data$InjectionUser & subject_data$RatingType == 'craving', 'Value']
  noninjection_craving_ratings <- subject_data[(subject_data$ImageSet == ImageSet) & (!subject_data$InjectionUser) & subject_data$RatingType == 'craving', 'Value']
  this_test <- t.test(injection_craving_ratings, noninjection_craving_ratings)
  merged_summary_values_sets[merged_summary_values_sets$Group == ImageSet, 'p_cravingdelta_injectionVnon'] <- this_test$p.value
}

library(dplyr)
library(data.table)

merged_summary_values <- rbindlist(list(merged_summary_values, merged_summary_values_sets), fill = TRUE, use.names = TRUE)
merged_summary_values <- data.frame(merged_summary_values)

###HSV code
#variables to write out, in a reasonable order
vars_to_save <- c("Group", "ImageSet", "File", "Category",
                  "allsubjects_valence_mean", "allsubjects_valence_sd", "allsubjects_arousal_mean", "allsubjects_arousal_sd", 
                  "allsubjects_craving_mean", "allsubjects_craving_sd", "allsubjects_typicality_mean", "allsubjects_typicality_sd",
                  "allsubjects_FractionNeither", "allsubjects_FractionMeth", "allsubjects_FractionOpioid", "allsubjects_FractionBoth", "allsubjects_MethToOpioid",
                  "injectionusers_valence_mean", "injectionusers_valence_sd", "injectionusers_arousal_mean", "injectionusers_arousal_sd", 
                  "injectionusers_craving_mean", "injectionusers_craving_sd", "injectionusers_typicality_mean", "injectionusers_typicality_sd",
                  "injectionusers_FractionNeither", "injectionusers_FractionMeth", "injectionusers_FractionOpioid", "injectionusers_FractionBoth", "injectionusers_MethToOpioid",
                  "noninjectionusers_valence_mean", "noninjectionusers_valence_sd", "noninjectionusers_arousal_mean", "noninjectionusers_arousal_sd", 
                  "noninjectionusers_craving_mean", "noninjectionusers_craving_sd", "noninjectionusers_typicality_mean", "noninjectionusers_typicality_sd",
                  "noninjectionusers_FractionNeither", "noninjectionusers_FractionMeth", "noninjectionusers_FractionOpioid", "noninjectionusers_FractionBoth",
                  "noninjectionusers_MethToOpioid",
                  "hue_mean", "hue_sd", "saturation_mean", "saturation_sd", "value_mean", "value_sd", "p_cravingdelta_injectionVnon")


write.csv(format(merged_summary_values[, vars_to_save], digits = 3), 'DCR_summaries_7-3-2019.csv', row.names = FALSE)

#only want image summaries, not category means, for example
plot_vals <- merged_summary_values[!is.na(merged_summary_values$File),]
#see how distance from neutral compares with mean valence
ggplot(data = plot_vals, aes(x = allsubjects_valence_mean, y = allsubjects_valence_fromneutral_mean, color = ImageSet)) + geom_point() + 
  xlab('valence') + ylab('valence distance to mean')+ geom_smooth(method = 'lm')  +
    scale_color_manual(values = cols)

print("#####Test for control images")
## [1] "#####Test for control images"
cor.test(plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == 'control'],plot_vals$allsubjects_valence_fromneutral_mean[plot_vals$ImageSet == 'control'] )
## 
##  Pearson's product-moment correlation
## 
## data:  plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == "control"] and plot_vals$allsubjects_valence_fromneutral_mean[plot_vals$ImageSet == plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == "control"] and     "control"]
## t = -0.0043273, df = 118, p-value = 0.9966
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1796269  0.1788557
## sample estimates:
##           cor 
## -0.0003983601
print("#####Test for meth images")
## [1] "#####Test for meth images"
cor.test(plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == 'meth'],plot_vals$allsubjects_valence_fromneutral_mean[plot_vals$ImageSet == 'meth'] )
## 
##  Pearson's product-moment correlation
## 
## data:  plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == "meth"] and plot_vals$allsubjects_valence_fromneutral_mean[plot_vals$ImageSet == plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == "meth"] and     "meth"]
## t = -6.3426, df = 118, p-value = 4.322e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6268160 -0.3572733
## sample estimates:
##       cor 
## -0.504225
print("#####Test for opioid images")
## [1] "#####Test for opioid images"
cor.test(plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == 'opioid'],plot_vals$allsubjects_valence_fromneutral_mean[plot_vals$ImageSet == 'opioid'] )
## 
##  Pearson's product-moment correlation
## 
## data:  plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == "opioid"] and plot_vals$allsubjects_valence_fromneutral_mean[plot_vals$ImageSet == plot_vals$allsubjects_valence_mean[plot_vals$ImageSet == "opioid"] and     "opioid"]
## t = -11.6, df = 118, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8039688 -0.6335581
## sample estimates:
##        cor 
## -0.7299109
scatter_data <- subject_data[subject_data$RatingType == c('valence'), c('id', 'ImageSet', 'ImageSetFile', 'RatingType', 'Value')]

names(scatter_data)[names(scatter_data) == 'Value'] <- 'Valence'
scatter_data$ValenceFromMean <- abs(scatter_data$Valence - 5)

#c <- cor(craving_data$Craving, typicality_data$Typicality)


ggplot(scatter_data, aes(x = Valence, y = ValenceFromMean, color = ImageSet)) + geom_point() + geom_smooth(method = 'lm')  +
    scale_color_manual(values = cols)

ggplot(scatter_data, aes(x = Valence, y = ValenceFromMean, color = ImageSet)) + geom_hex() + geom_smooth(method = 'lm')  +
    scale_color_manual(values = cols)
## Warning: package 'hexbin' was built under R version 3.5.3

Scatter plot of individual craving/typicality ratings, just within drug cues

scatter_data <- subject_data[subject_data$ImageSet %in% c('meth', 'opioid'), c('id', 'ImageSet', 'ImageSetFile', 'RatingType', 'Value')]

craving_data <- scatter_data[scatter_data$RatingType == 'craving',]
typicality_data <- scatter_data[scatter_data$RatingType == 'typicality',]

craving_data$Craving <- craving_data$Value
typicality_data$Typicality <- typicality_data$Value

correlation <- cor(craving_data$Craving, typicality_data$Typicality)

scatter_data <- merge(craving_data[, c('id', 'ImageSet', 'ImageSetFile', 'Craving')], typicality_data[, c('id', 'ImageSet', 'ImageSetFile', 'Typicality')])

ggplot(scatter_data, aes(x = Craving, y = Typicality, color = ImageSet)) + geom_point() + geom_smooth(method = 'lm') + 
  ggtitle(paste0('r = ', round(correlation, 2), '\nBut not appropriate, since multiple measures per subject')) +
    scale_color_manual(values = cols)

mean_scatter_data <- merged_summary_values[merged_summary_values$ImageSet %in% c('meth', 'opioid'),]
correlation <- cor(mean_scatter_data$allsubjects_craving_mean, mean_scatter_data$allsubjects_typicality_mean)
ggplot(mean_scatter_data, aes(x = allsubjects_craving_mean, y = allsubjects_typicality_mean, color = ImageSet)) + geom_point() + geom_smooth(method = 'lm')  +
  scale_y_continuous(breaks = c(0, 25, 50, 75, 100), limits = c(0, 100)) +
  scale_x_continuous(breaks = c(0, 25, 50, 75, 100), limits = c(0, 100)) +
  xlab('Craving Mean') + ylab ('Typicality Mean') + ggtitle(paste0('r = ', round(correlation, 2))) +
  scale_color_manual(values = cols)

Histograms for each image set

make_plotset <- function(data_to_plot, main_title, filename, group){
  p1 <- ggplot(data_to_plot, aes_string(x = 'allsubjects_valence_mean', fill = 'ImageSet')) + 
    geom_histogram(binwidth = 0.5, alpha = 0.5, position = 'identity') + theme(legend.position="none") +
    scale_color_manual(values = cols)
  p2 <- ggplot(data_to_plot, aes_string(x = 'allsubjects_arousal_mean', fill = 'ImageSet')) + 
    geom_histogram(binwidth = 0.5, alpha = 0.5, position = 'identity') + theme(legend.position="none") +
    scale_color_manual(values = cols)
  p3 <- ggplot(data_to_plot, aes_string(x = 'allsubjects_craving_mean', fill = 'ImageSet')) + 
    geom_histogram(binwidth = 10, alpha = 0.5, position = 'identity') + theme(legend.position="none") +
    scale_color_manual(values = cols)
  p4 <- ggplot(data_to_plot, aes_string(x = 'allsubjects_typicality_mean', fill = 'ImageSet')) + 
    geom_histogram(binwidth = 10, alpha = 0.5, position = 'identity') + theme(legend.position="none") +
    scale_color_manual(values = cols)
  p5 <- ggplot(data_to_plot, aes_string(x = 'hue_mean', fill = 'ImageSet')) + 
    geom_histogram(binwidth = 0.1, alpha = 0.5, position = 'identity') + theme(legend.position="none") +
    scale_color_manual(values = cols)
  p6 <- ggplot(data_to_plot, aes_string(x = 'saturation_mean', fill = 'ImageSet')) + 
    geom_histogram(binwidth = 0.1, alpha = 0.5, position = 'identity') + theme(legend.position="none") +
    scale_color_manual(values = cols)
  p7 <- ggplot(data_to_plot, aes_string(x = 'value_mean', fill = 'ImageSet')) + 
    geom_histogram(binwidth = 0.1, alpha = 0.5, position = 'identity') +
    scale_color_manual(values = cols)
  #png(filename, width = 1400, height = 400)
  #grid.arrange(p1, p2, p3, p4, p5, p6, p7, ncol = 7, top = textGrob(main_title, gp=gpar(fontsize = 50)))
  #dev.off()
  print(p1)
  print(p2)
  print(p3)
  print(p4)
  print(p5)
  print(p6)
  print(p7)
}
library(sm)
## Warning: package 'sm' was built under R version 3.5.3
## Package 'sm', version 2.2-5.6: type help(sm) for summary information
one_density_plot <- function(var_to_plot, label){

  sm.density.compare(merged_summary_values[, var_to_plot], as.factor(merged_summary_values$ImageSet), lwd = 5)
  title(main=label)
  colfill<-c(2:(2+length(levels(as.factor(merged_summary_values$ImageSet))))) 
  legend('topright', levels(as.factor(merged_summary_values$ImageSet)), fill=colfill)
}



source('R_rainclouds.R')
library(cowplot)
## Warning: package 'cowplot' was built under R version 3.5.3
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
## 
##     get_legend
## The following object is masked from 'package:ggplot2':
## 
##     ggsave
one_raincloud_plot <- function(var_to_plot, label){


  p6 <- ggplot(subject_data[subject_data$RatingType == 'arousal',],aes(x=ImageSet,y=Value, fill = ImageSet, colour = ImageSet))+
    geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim =  FALSE)+
    #geom_point(position = position_jitter(width = .15), size = .25)+
    geom_point(position = position_jitter(width = 1), size = .25)+
    geom_boxplot(aes(x = as.numeric(ImageSet)+0.25, y = Value),outlier.shape = NA,
    alpha = 0.3, width = .1, colour = "BLACK") +
    ylab('Arousal')+xlab('ImageSet')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
    colour = FALSE) +
    ggtitle("Figure R6: Change in Colour Palette") +
    scale_color_manual(values = cols) +
    scale_fill_manual(values = cols)
    #scale_colour_brewer(palette = "Dark2")+
    #scale_fill_brewer(palette = "Dark2")+
  ggsave(paste0(label, '.png'), width = w, height = h)
  p6
}

subject_data$ImageSetCap <- as.character(subject_data$ImageSet)
subject_data$ImageSetCap[subject_data$ImageSetCap == 'meth'] <- 'Meth'
subject_data$ImageSetCap[subject_data$ImageSetCap == 'opioid'] <- 'Opioid'
subject_data$ImageSetCap[subject_data$ImageSetCap == 'control'] <- 'Control'
subject_data$ImageSetCap <- factor(subject_data$ImageSetCap, levels = c('Opioid', 'Meth', 'Control'))

p1 <- ggplot(subject_data[subject_data$RatingType == 'craving',],aes(x=ImageSetCap,y=Value, fill = ImageSetCap, colour = ImageSetCap))+
  geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim =  FALSE)+
  geom_point(position = position_jitter(width = .15), size = .25)+
  geom_boxplot(aes(x = as.numeric(ImageSetCap)+0.25, y = Value),outlier.shape = NA,
  alpha = 0.3, width = .1, colour = "BLACK") +
  ylab('')+xlab('')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
  colour = FALSE)  +
  scale_color_manual(values = cols) +
  scale_fill_manual(values = cols) +
  ggtitle("Craving") + scale_y_continuous(breaks = c(0, 25, 50, 75, 100), limits = c(0, 100))

p2 <- ggplot(subject_data[subject_data$RatingType == 'arousal',],aes(x=ImageSetCap,y=Value, fill = ImageSetCap, colour = ImageSetCap))+
  geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim =  FALSE)+
  geom_point(position = position_jitter(width = .15), size = .25)+
  geom_boxplot(aes(x = as.numeric(ImageSetCap)+0.25, y = Value),outlier.shape = NA,
  alpha = 0.3, width = .1, colour = "BLACK") +
  ylab('')+xlab('')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
  colour = FALSE)  +
  scale_color_manual(values = cols) +
  scale_fill_manual(values = cols) +
  ggtitle("Arousal") + scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9))

p3 <- ggplot(subject_data[subject_data$RatingType == 'valence',],aes(x=ImageSetCap,y=Value, fill = ImageSetCap, colour = ImageSetCap))+
  geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim =  FALSE)+
  geom_point(position = position_jitter(width = .15), size = .25)+
  geom_boxplot(aes(x = as.numeric(ImageSetCap)+0.25, y = Value),outlier.shape = NA,
  alpha = 0.3, width = .1, colour = "BLACK") +
  ylab('')+xlab('')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
  colour = FALSE) +
  scale_color_manual(values = cols) +
  scale_fill_manual(values = cols) +
  ggtitle("Valence") + scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9))



png('rainclouds.png', width = 1400, height = 400)
print(annotate_figure(ggarrange(p1, p2, p3, ncol = 3, nrow = 1, common.legend = TRUE, legend = 'right'), top = text_grob('Rainclouds', color = 'red')))
## Warning: Removed 168 rows containing missing values (geom_flat_violin).
## Warning: Removed 1858 rows containing missing values (geom_point).
## Warning: Removed 380 rows containing missing values (geom_flat_violin).
## Warning: Removed 2096 rows containing missing values (geom_point).
## Warning: Removed 252 rows containing missing values (geom_flat_violin).
## Warning: Removed 464 rows containing missing values (geom_point).
## Warning: Removed 168 rows containing missing values (geom_flat_violin).
## Warning: Removed 1858 rows containing missing values (geom_point).
## Warning: Removed 380 rows containing missing values (geom_flat_violin).
## Warning: Removed 2096 rows containing missing values (geom_point).
## Warning: Removed 252 rows containing missing values (geom_flat_violin).
## Warning: Removed 464 rows containing missing values (geom_point).
dev.off()
## png 
##   2
p4 <- ggplot(subject_data,aes(x=ImageSetCap,y=millisecondsOnPage / 1000, fill = ImageSetCap, colour = ImageSetCap))+
  geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim =  FALSE)+
  geom_point(position = position_jitter(width = .15), size = .25)+
  geom_boxplot(aes(x = as.numeric(ImageSetCap)+0.25, y = millisecondsOnPage / 1000),outlier.shape = NA,
  alpha = 0.3, width = .1, colour = "BLACK") +
  ylab('')+xlab('')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
  colour = FALSE) +
  scale_color_manual(values = cols) +
  scale_fill_manual(values = cols) +
  ggtitle("SecondsOnPage") #+ scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9))

p5 <- ggplot(subject_data,aes(x=ImageSetCap,y=responseDelay / 1000, fill = ImageSetCap, colour = ImageSetCap))+
  geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim =  FALSE)+
  geom_point(position = position_jitter(width = .15), size = .25)+
  geom_boxplot(aes(x = as.numeric(ImageSetCap)+0.25, y = responseDelay / 1000),outlier.shape = NA,
  alpha = 0.3, width = .1, colour = "BLACK") +
  ylab('')+xlab('')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
  colour = FALSE) +
  scale_color_manual(values = cols) +
  scale_fill_manual(values = cols) +
  ggtitle("ResponseDelay") #+ scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9))

png('RTrainclouds.png', width = 1400, height = 400)
print(annotate_figure(ggarrange(p4, p5, ncol = 1, nrow = 2, common.legend = TRUE, legend = 'right'), top = text_grob('RTRainclouds', color = 'red')))
dev.off()
## png 
##   2
p1 <- ggplot(subject_data[subject_data$RatingType == 'craving',], 
            aes(x = Value, y = millisecondsOnPage / 1000,color = ImageSet)) + geom_point() +
  xlab('Craving') + ylab('ScondsOnPage') + #scale_x_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) + 
    #scale_y_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) + 
  theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
  geom_smooth(method = 'lm') +
  scale_color_manual(values = cols) +
  scale_fill_manual(values = cols)
p2 <- ggplot(subject_data[subject_data$RatingType == 'craving',], 
            aes(x = Value, y = responseDelay / 1000,color = ImageSet)) + geom_point() +
  xlab('Craving') + ylab('ResponseDelay') + #scale_x_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) + 
    #scale_y_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) + 
  theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
  geom_smooth(method = 'lm') +
  scale_color_manual(values = cols) +
  scale_fill_manual(values = cols)

png('RTCravingScatter.png', width = 400, height = 800)
print(annotate_figure(ggarrange(p1, p2, ncol = 1, nrow = 2, common.legend = TRUE, legend = 'right'), top = text_grob('RTScatters', color = 'red')))
dev.off()
## png 
##   2
craving_data <- subject_data[subject_data$RatingType == 'craving',]

#set threshold of 20 seconds, remove 0.6% of data
sum(craving_data$millisecondsOnPage > 10000) / nrow(craving_data)
## [1] 0.02434343
one_density_plot('allsubjects_valence_mean', 'Valence')

one_density_plot('allsubjects_arousal_mean', 'Arousal')

one_density_plot('allsubjects_craving_mean', 'Craving')

one_density_plot('allsubjects_typicality_mean', 'Typicality')

one_density_plot('hue_mean', 'Hue')

one_density_plot('saturation_mean', 'Saturation')

one_density_plot('value_mean', 'Value')

Correlations Between Mean Ratings for Each Image Set

library(corrplot)
## corrplot 0.84 loaded
library(PerformanceAnalytics)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:data.table':
## 
##     first, last
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
cols_to_use <-  c('allsubjects_valence_mean', 'allsubjects_valence_fromneutral_mean', 'allsubjects_arousal_mean', 
                  'allsubjects_craving_mean')

one_corrplot <- function(dset, title){
  
  col2 <- colorRampPalette(rev(c("#67001F", "#B2182B", "#D6604D", "#F4A582",
                           "#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
                           "#4393C3", "#2166AC", "#053061")))

  
  corvals <- cor(dset)
  rownames(corvals) <- colnames(corvals) <- c('valence', 'abs_valence', 'arousal', 'craving')
  corrplot.mixed(corvals, upper = 'ellipse', lower = 'number', title = title, mar=c(0,0,2,0), 
                 upper.col = col2(50), lower.col = col2(50))
  chart.Correlation(dset)
  #returns z-transformed correlations, for easier comparisons between matrices
  return(atanh(corvals))
}

p <- ggplot(merged_summary_values[merged_summary_values$ImageSet %in% c('meth', 'opioid'),], aes(x = allsubjects_craving_mean, y = allsubjects_typicality_mean,
                                                                                           color = ImageSet)) + geom_point() +
  xlab('Craving') + ylab('Typicality') + scale_x_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) + 
    scale_y_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) + 
  theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
  geom_smooth(method = 'lm') +
  scale_color_manual(values = cols)



#wide_subject_data <- dcast(subject_data, formula = id ~ RatingType + ImageSetFile, value.var = 'Value')


p <- ggplot(subject_data[subject_data$ImageSet %in% c('meth', 'opioid'),], 
            aes(x = allsubjects_craving_mean, y = allsubjects_typicality_mean,color = ImageSet)) + geom_point() +
  xlab('Craving') + ylab('Typicality') + scale_x_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) + 
    scale_y_continuous(breaks = c(60, 70, 80, 90, 100), limits = c(60, 100)) + 
  theme(axis.text = element_text(size = 18), plot.title = element_text(size = 24), legend.text=element_text(size=14)) +
  geom_smooth(method = 'lm') +
  scale_color_manual(values = cols)





drug_zs <- one_corrplot(merged_summary_values[merged_summary_values$ImageSet %in% c('meth', 'opioid'), cols_to_use], 'Meth/Opioid')

meth_zs <- one_corrplot(merged_summary_values[merged_summary_values$ImageSet %in% c('meth'), cols_to_use], 'Meth')

opioid_zs <- one_corrplot(merged_summary_values[merged_summary_values$ImageSet %in% c('opioid'), cols_to_use], 'Opioid')

control_zs <- one_corrplot(merged_summary_values[merged_summary_values$ImageSet %in% c('control'), cols_to_use], 'control')

#to compare correlation z-scores--here we use N1=N2=120, the number of images
#Zobserved = (z1 - z2) / (square root of [ (1 / N1 - 3) + (1 / N2 - 3) ])
meth_opioid_deltazs <- (meth_zs - opioid_zs) / sqrt(1/(120-3) + 1/(120-3))
meth_opioid_deltaps <- 2*pnorm(-abs(meth_opioid_deltazs))

col2 <- colorRampPalette(rev(c("#67001F", "#B2182B", "#D6604D", "#F4A582",
                           "#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
                           "#4393C3", "#2166AC", "#053061")))


corrplot.mixed(meth_opioid_deltazs, upper = 'color', lower = 'number', title = 'Delta Z Score (Meth - Opioid)', mar=c(0,0,2,0), is.corr = FALSE,
               p.mat = meth_opioid_deltaps, sig.level = 1, lower.col = col2(50), upper.col = col2(50))

corrplot.mixed(meth_opioid_deltaps, upper = 'color', lower = 'number', title = 'Delta r p-values (Meth - Opioid)', mar=c(0,0,2,0), is.corr = FALSE,
               lower.col = col2(50), upper.col = col2(50))

#               p.mat = meth_opioid_deltaps, sig.level = 0.01, lower.col = col2(50), upper.col = col2(50))
print('Matrix of p-values for Correlation Differences')
## [1] "Matrix of p-values for Correlation Differences"
print(meth_opioid_deltaps)
##                 valence abs_valence     arousal   craving
## valence             NaN 0.003167665 0.005306939 0.0232624
## abs_valence 0.003167665         NaN 0.130490960 0.7358989
## arousal     0.005306939 0.130490960         NaN 0.2669365
## craving     0.023262397 0.735898874 0.266936547       NaN
print('Matrix of fdr corrected p-values')
## [1] "Matrix of fdr corrected p-values"
p.adjust(meth_opioid_deltaps[upper.tri(meth_opioid_deltaps)], method = 'fdr')
## [1] 0.01592082 0.01592082 0.19573644 0.04652479 0.73589887 0.32032386

Craving LMEs

meth_craving ~ Age + sex + age_of_meth_onset + duration_of_meth_use + cost_of_meth + days_of_meth + duration_of_current_abstinance

#first, corrplot of these variables
clinical_columns <- c('Age', 'SexMale', 'InjectionUser', 'MethOnsetAge', 'MethMonthsUsed', 'MethDollars',
                   'MethDaysLastMonth','MethDaysAbstinant', 'OpioidOnsetAge', 'OpioidMonthsUsed', 
                   'OpioidDaysLastMonth', 'OpioidDaysAbstinant', 'HeroinOnsetAge', 'HeroinMonthsUsed', 
                   'HeroinDaysLastMonth', 'HeroinDaysAbstinant', 'OpioidOrHeroinOnsetAge', 
                   'OpioidOrHeroinMonthsUsed', 'OpioidOrHeroinDollars', 'OpioidOrHeroinDaysLastMonth','OpioidOrHeroinDaysAbstinant') 
col2 <- colorRampPalette(rev(c("#67001F", "#B2182B", "#D6604D", "#F4A582",
                           "#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
                           "#4393C3", "#2166AC", "#053061")))

clinical_data$SexMale <- 0
clinical_data$SexMale[clinical_data$Sex == 'Male'] <- 1

corvals <- cor(clinical_data[, clinical_columns], use = 'pairwise.complete')
png('ClinicalCorrelations.png', width = 1000, height = 1000)
corrplot.mixed(corvals, upper = 'ellipse', lower = 'number', title = 'Clinical Variable Correlations', mar=c(0,0,2,0), 
               upper.col = col2(50), lower.col = col2(50), tl.pos = 'lt')
dev.off()
## png 
##   2
png('ClinicalScatterplotMatrix.png', width = 1000, height = 1000)
chart.Correlation(clinical_data[, clinical_columns])
dev.off()
## png 
##   2
library(tableone)
CreateTableOne(vars = clinical_columns, data = clinical_data)
##                                          
##                                           Overall          
##   n                                            28          
##   Age (mean (sd))                           37.71 (8.11)   
##   SexMale (mean (sd))                        0.57 (0.50)   
##   InjectionUser = TRUE (%)                     21 (75.0)   
##   MethOnsetAge (mean (sd))                  19.42 (5.40)   
##   MethMonthsUsed (mean (sd))               117.96 (96.28)  
##   MethDollars (mean (sd))                 2928.93 (2540.84)
##   MethDaysLastMonth (mean (sd))             21.00 (12.29)  
##   MethDaysAbstinant (mean (sd))            729.32 (858.29) 
##   OpioidOnsetAge (mean (sd))                20.82 (6.39)   
##   OpioidMonthsUsed (mean (sd))             115.32 (94.96)  
##   OpioidDaysLastMonth (mean (sd))           16.71 (12.91)  
##   OpioidDaysAbstinant (mean (sd))          751.04 (664.18) 
##   HeroinOnsetAge (mean (sd))                26.79 (7.73)   
##   HeroinMonthsUsed (mean (sd))              31.43 (71.31)  
##   HeroinDaysLastMonth (mean (sd))           10.61 (13.42)  
##   HeroinDaysAbstinant (mean (sd))         1492.20 (2022.00)
##   OpioidOrHeroinOnsetAge (mean (sd))        20.54 (6.29)   
##   OpioidOrHeroinMonthsUsed (mean (sd))     117.89 (92.67)  
##   OpioidOrHeroinDollars (mean (sd))       1773.21 (1711.63)
##   OpioidOrHeroinDaysLastMonth (mean (sd))   22.36 (11.18)  
##   OpioidOrHeroinDaysAbstinant (mean (sd))  914.43 (1332.91)
library(nlme)
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
library(sjmisc)
## Warning: package 'sjmisc' was built under R version 3.5.3
library(sjPlot)
## Warning: package 'sjPlot' was built under R version 3.5.3
## 
## Attaching package: 'sjPlot'
## The following objects are masked from 'package:cowplot':
## 
##     plot_grid, save_plot
library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
## 
##     lmList
library(lmerTest)
## Warning: package 'lmerTest' was built under R version 3.5.3
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
this_lme <- lmer(Value ~ Age + Sex + MethOnsetAge + MethMonthsUsed + MethDollars + MethDaysLastMonth + MethDaysAbstinant + Time * Visit +(1 | id), 
              data = subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet == 'meth',])
## Warning: Some predictor variables are on very different scales: consider rescaling
## Warning: Some predictor variables are on very different scales: consider rescaling
print(plot_model(this_lme, title = 'meth craving', show.values = TRUE, type = 'std', order.terms = c(1,2,8,9,10,3,4,6,5,7)))

print(summary(this_lme))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Age + Sex + MethOnsetAge + MethMonthsUsed + MethDollars +  
##     MethDaysLastMonth + MethDaysAbstinant + Time * Visit + (1 |      id)
##    Data: subject_data[subject_data$RatingType == "craving" & subject_data$ImageSet ==  
##     "meth", ]
## 
## REML criterion at convergence: 26081.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.9873 -0.1122  0.1452  0.5159  2.9383 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 132.9    11.53   
##  Residual             204.4    14.30   
## Number of obs: 3180, groups:  id, 27
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        8.467e+01  1.915e+01  1.907e+01   4.421 0.000291 ***
## Age                4.888e-01  3.550e-01  1.901e+01   1.377 0.184622    
## SexMale           -8.208e+00  5.911e+00  1.901e+01  -1.389 0.181023    
## MethOnsetAge       7.246e-02  5.766e-01  1.902e+01   0.126 0.901319    
## MethMonthsUsed     9.298e-03  3.612e-02  1.900e+01   0.257 0.799619    
## MethDollars        1.451e-03  1.191e-03  1.904e+01   1.218 0.238067    
## MethDaysLastMonth -5.150e-01  2.326e-01  1.902e+01  -2.213 0.039282 *  
## MethDaysAbstinant -4.520e-03  8.542e-03  1.900e+01  -0.529 0.602841    
## Time              -2.011e-02  6.679e-03  3.150e+03  -3.011 0.002625 ** 
## VisitV2            3.247e-01  1.029e+00  3.151e+03   0.316 0.752341    
## Time:VisitV2      -4.525e-02  9.591e-03  3.150e+03  -4.717 2.49e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Age    SexMal MthOnA MthMnU MthDll MthDLM MthDyA Time   VistV2
## Age         -0.526                                                               
## SexMale     -0.061 -0.002                                                        
## MethOnsetAg -0.570 -0.154 -0.386                                                 
## MthMnthsUsd -0.222 -0.384 -0.330  0.580                                          
## MethDollars -0.337  0.114 -0.009  0.166 -0.053                                   
## MthDysLstMn -0.155 -0.110  0.346 -0.016 -0.074 -0.433                            
## MthDysAbstn -0.507  0.072  0.334  0.051  0.101  0.386  0.058                     
## Time        -0.032  0.000 -0.002  0.002  0.002 -0.002 -0.001 -0.004              
## VisitV2     -0.022 -0.004  0.002 -0.002  0.003 -0.003  0.001 -0.003  0.602       
## Time:VistV2  0.021  0.002  0.000  0.000 -0.002  0.000  0.002  0.002 -0.695 -0.868
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
print(anova(this_lme))
## Type III Analysis of Variance Table with Satterthwaite's method
##                    Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## Age                 387.3   387.3     1   19.01  1.8951   0.18462    
## Sex                 394.0   394.0     1   19.01  1.9281   0.18102    
## MethOnsetAge          3.2     3.2     1   19.02  0.0158   0.90132    
## MethMonthsUsed       13.5    13.5     1   19.00  0.0663   0.79962    
## MethDollars         303.2   303.2     1   19.04  1.4838   0.23807    
## MethDaysLastMonth  1001.3  1001.3     1   19.02  4.8996   0.03928 *  
## MethDaysAbstinant    57.2    57.2     1   19.00  0.2800   0.60284    
## Time              16175.7 16175.7     1 3150.40 79.1509 < 2.2e-16 ***
## Visit                20.4    20.4     1 3150.88  0.0996   0.75234    
## Time:Visit         4547.9  4547.9     1 3150.24 22.2536 2.493e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
scaled_data <- subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet == 'meth',]
scaled_data[, c('Age', 'MethOnsetAge',  'MethMonthsUsed', 'MethDollars',  'MethDaysLastMonth',  'MethDaysAbstinant')] <-
  scale(scaled_data[, c('Age', 'MethOnsetAge',  'MethMonthsUsed', 'MethDollars',  'MethDaysLastMonth',  'MethDaysAbstinant')])
this_lme <- lmer(Value ~ Age + Sex + MethOnsetAge + MethMonthsUsed + MethDollars + MethDaysLastMonth + MethDaysAbstinant + Time * Visit + (1 | id), 
              data = scaled_data)
print(plot_model(this_lme, title = 'meth craving scaled predictors', show.values = TRUE, type = 'std', order.terms = c(1,2,8,9,10,3,4,6,5,7)))

print(summary(this_lme))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Age + Sex + MethOnsetAge + MethMonthsUsed + MethDollars +  
##     MethDaysLastMonth + MethDaysAbstinant + Time * Visit + (1 |      id)
##    Data: scaled_data
## 
## REML criterion at convergence: 26030.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.9873 -0.1122  0.1452  0.5159  2.9383 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 132.9    11.53   
##  Residual             204.4    14.30   
## Number of obs: 3180, groups:  id, 27
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        9.582e+01  4.158e+00  2.003e+01  23.041 6.85e-16 ***
## Age                3.926e+00  2.852e+00  1.901e+01   1.377  0.18462    
## SexMale           -8.208e+00  5.911e+00  1.901e+01  -1.389  0.18102    
## MethOnsetAge       3.795e-01  3.020e+00  1.902e+01   0.126  0.90132    
## MethMonthsUsed     8.708e-01  3.383e+00  1.900e+01   0.257  0.79962    
## MethDollars        3.651e+00  2.998e+00  1.904e+01   1.218  0.23807    
## MethDaysLastMonth -6.240e+00  2.819e+00  1.902e+01  -2.213  0.03928 *  
## MethDaysAbstinant -3.834e+00  7.246e+00  1.900e+01  -0.529  0.60284    
## Time              -2.011e-02  6.679e-03  3.150e+03  -3.011  0.00262 ** 
## VisitV2            3.247e-01  1.029e+00  3.151e+03   0.316  0.75234    
## Time:VisitV2      -4.525e-02  9.591e-03  3.150e+03  -4.717 2.49e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Age    SexMal MthOnA MthMnU MthDll MthDLM MthDyA Time   VistV2
## Age          0.076                                                               
## SexMale     -0.761 -0.002                                                        
## MethOnsetAg  0.371 -0.154 -0.386                                                 
## MthMnthsUsd  0.336 -0.384 -0.330  0.580                                          
## MethDollars  0.136  0.114 -0.009  0.166 -0.053                                   
## MthDysLstMn -0.298 -0.110  0.346 -0.016 -0.074 -0.433                            
## MthDysAbstn  0.047  0.072  0.334  0.051  0.101  0.386  0.058                     
## Time        -0.148  0.000 -0.002  0.002  0.002 -0.002 -0.001 -0.004              
## VisitV2     -0.122 -0.004  0.002 -0.002  0.003 -0.003  0.001 -0.003  0.602       
## Time:VistV2  0.104  0.002  0.000  0.000 -0.002  0.000  0.002  0.002 -0.695 -0.868
print(anova(this_lme))
## Type III Analysis of Variance Table with Satterthwaite's method
##                    Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## Age                 387.3   387.3     1   19.01  1.8951   0.18462    
## Sex                 394.0   394.0     1   19.01  1.9281   0.18102    
## MethOnsetAge          3.2     3.2     1   19.02  0.0158   0.90132    
## MethMonthsUsed       13.5    13.5     1   19.00  0.0663   0.79962    
## MethDollars         303.2   303.2     1   19.04  1.4838   0.23807    
## MethDaysLastMonth  1001.3  1001.3     1   19.02  4.8996   0.03928 *  
## MethDaysAbstinant    57.2    57.2     1   19.00  0.2800   0.60284    
## Time              16175.7 16175.7     1 3150.40 79.1509 < 2.2e-16 ***
## Visit                20.4    20.4     1 3150.88  0.0996   0.75234    
## Time:Visit         4547.9  4547.9     1 3150.24 22.2536 2.493e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
this_lme <- lmer(Value ~ Age + Sex + OpioidOrHeroinOnsetAge + OpioidOrHeroinMonthsUsed + OpioidOrHeroinDollars + 
                   OpioidOrHeroinDaysLastMonth + OpioidOrHeroinDaysAbstinant + Time * Visit + (1 | id), 
              data = subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet == 'opioid',])
## Warning: Some predictor variables are on very different scales: consider rescaling

## Warning: Some predictor variables are on very different scales: consider rescaling
print(plot_model(this_lme, title = 'opioid craving', show.values = TRUE, type = 'std', order.terms = c(1,2,8, 9, 10, 3,4,6,5,7)))

print(summary(this_lme))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Age + Sex + OpioidOrHeroinOnsetAge + OpioidOrHeroinMonthsUsed +  
##     OpioidOrHeroinDollars + OpioidOrHeroinDaysLastMonth + OpioidOrHeroinDaysAbstinant +  
##     Time * Visit + (1 | id)
##    Data: subject_data[subject_data$RatingType == "craving" & subject_data$ImageSet ==  
##     "opioid", ]
## 
## REML criterion at convergence: 25942.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.3037 -0.1131  0.1809  0.4864  3.1743 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 139.1    11.80   
##  Residual             195.4    13.98   
## Number of obs: 3180, groups:  id, 27
## 
## Fixed effects:
##                               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                  8.189e+01  1.626e+01  1.905e+01   5.036 7.28e-05 ***
## Age                          2.123e-01  4.675e-01  1.900e+01   0.454 0.654880    
## SexMale                     -5.248e+00  5.059e+00  1.903e+01  -1.037 0.312565    
## OpioidOrHeroinOnsetAge       1.731e-01  5.309e-01  1.904e+01   0.326 0.747881    
## OpioidOrHeroinMonthsUsed     5.140e-02  3.799e-02  1.902e+01   1.353 0.191893    
## OpioidOrHeroinDollars        1.204e-03  1.784e-03  1.904e+01   0.675 0.507923    
## OpioidOrHeroinDaysLastMonth -2.161e-01  2.839e-01  1.900e+01  -0.761 0.455866    
## OpioidOrHeroinDaysAbstinant  3.980e-04  2.047e-03  1.901e+01   0.194 0.847885    
## Time                        -3.784e-02  6.548e-03  3.150e+03  -5.779 8.27e-09 ***
## VisitV2                     -3.337e+00  9.908e-01  3.151e+03  -3.368 0.000767 ***
## Time:VisitV2                -1.299e-02  9.328e-03  3.150e+03  -1.392 0.164007    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Age    SexMal OpOHOA OpOHMU OpdOHD OOHDLM OpOHDA Time   VistV2
## Age         -0.756                                                               
## SexMale     -0.067  0.019                                                        
## OpdOrHrnOnA  0.085 -0.616 -0.142                                                 
## OpdOrHrnMnU  0.128 -0.396 -0.332  0.405                                          
## OpdOrHrnDll  0.057  0.117  0.209 -0.417 -0.245                                   
## OpdOrHrnDLM -0.633  0.487  0.057 -0.264 -0.393 -0.157                            
## OpdOrHrnDyA -0.177 -0.044 -0.157 -0.007  0.225  0.020  0.235                     
## Time        -0.034 -0.001  0.001 -0.002 -0.001  0.001 -0.002  0.000              
## VisitV2     -0.026 -0.001  0.002 -0.003 -0.003  0.004 -0.002 -0.001  0.603       
## Time:VistV2  0.023  0.003  0.001 -0.001  0.001 -0.001  0.002 -0.001 -0.702 -0.863
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
print(anova(this_lme))
## Type III Analysis of Variance Table with Satterthwaite's method
##                              Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## Age                            40.3    40.3     1   19.00  0.2062 0.6548801    
## Sex                           210.3   210.3     1   19.03  1.0761 0.3125645    
## OpioidOrHeroinOnsetAge         20.8    20.8     1   19.04  0.1064 0.7478808    
## OpioidOrHeroinMonthsUsed      357.7   357.7     1   19.02  1.8309 0.1918931    
## OpioidOrHeroinDollars          89.0    89.0     1   19.04  0.4553 0.5079235    
## OpioidOrHeroinDaysLastMonth   113.2   113.2     1   19.00  0.5795 0.4558659    
## OpioidOrHeroinDaysAbstinant     7.4     7.4     1   19.01  0.0378 0.8478851    
## Time                        17646.5 17646.5     1 3150.25 90.3210 < 2.2e-16 ***
## Visit                        2215.8  2215.8     1 3150.89 11.3411 0.0007672 ***
## Time:Visit                    378.6   378.6     1 3150.27  1.9378 0.1640067    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

LMEs Of Image Ratings

We will look drug and neutral images separately. For drug images, the model will be:

Value ~ Time * Visit + (Age * Sex * ImageSet)^2 + (1 | id)

And for neutral images it will be:

Value ~ Time * Visit + (Age * Sex) + (1 | id)

for valence, arousal, and craving separately

summarize_one_lme_drug <- function(dset, RatingType){
  this_lme <- lmer(Value ~ Time * Visit + (Age * Sex * ImageSet)^2 + (1 | id), 
                data = dset[dset$RatingType == RatingType,])
  print(plot_model(this_lme, title = RatingType, show.values = TRUE, type = 'std'))
  print(summary(this_lme))
  #for using nlme instead of lme4--but can't get standardized betas in the forest plot
  #that_lme <- lme(fixed = Value ~ Time * Visit + Age + Sex + ImageSet, random = ~ 1 | id, 
  #              data = dset[dset$RatingType == RatingType,])
  #print(summary(that_lme))
  print(anova(this_lme))
}
summarize_one_lme_neutral <- function(dset, RatingType){
  this_lme <- lmer(Value ~ Time * Visit + (Age * Sex) + (1 | id), 
                data = dset[dset$RatingType == RatingType,])
  print(plot_model(this_lme, title = RatingType, show.values = TRUE, type = 'std'))
  print(summary(this_lme))
  print(anova(this_lme))
}

Drug Cues

summarize_one_lme_drug(subject_data[subject_data$ImageSet %in% c('meth', 'opioid'),], 'valence')

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex * ImageSet)^2 + (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 18920.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.8998 -0.3812  0.0481  0.4229  6.6748 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 2.041    1.4287  
##  Residual             0.995    0.9975  
## Number of obs: 6597, groups:  id, 28
## 
## Fixed effects:
##                              Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                 3.863e+00  1.687e+00  2.412e+01   2.290   0.0311 *  
## Time                        2.065e-03  3.230e-04  6.562e+03   6.393 1.74e-10 ***
## VisitV2                     3.917e-01  4.943e-02  6.562e+03   7.925 2.66e-15 ***
## Age                        -3.369e-03  4.350e-02  2.410e+01  -0.077   0.9389    
## SexMale                    -5.422e+00  2.682e+00  2.410e+01  -2.022   0.0545 .  
## ImageSetopioid             -7.434e-02  1.518e-01  6.562e+03  -0.490   0.6243    
## Time:VisitV2               -4.095e-03  4.632e-04  6.562e+03  -8.842  < 2e-16 ***
## Age:SexMale                 1.661e-01  6.960e-02  2.410e+01   2.386   0.0252 *  
## Age:ImageSetopioid          2.162e-03  3.915e-03  6.562e+03   0.552   0.5809    
## SexMale:ImageSetopioid      5.365e-01  2.421e-01  6.562e+03   2.216   0.0267 *  
## Age:SexMale:ImageSetopioid -1.445e-02  6.273e-03  6.562e+03  -2.303   0.0213 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Time   VistV2 Age    SexMal ImgStp Tm:VV2 Ag:SxM Ag:ImS SxM:IS
## Time        -0.018                                                               
## VisitV2     -0.013  0.602                                                        
## Age         -0.969  0.000 -0.001                                                 
## SexMale     -0.629  0.000  0.000  0.610                                          
## ImageSetopd -0.045  0.001 -0.013  0.044  0.028                                   
## Time:VistV2  0.011 -0.697 -0.865  0.001  0.001  0.015                            
## Age:SexMale  0.606  0.000  0.000 -0.625 -0.979 -0.027 -0.001                     
## Age:ImgStpd  0.044  0.002  0.018 -0.045 -0.027 -0.969 -0.021  0.028              
## SxMl:ImgStp  0.028  0.005  0.012 -0.027 -0.045 -0.627 -0.013  0.044  0.608       
## Ag:SxMl:ImS -0.027 -0.006 -0.016  0.028  0.044  0.605  0.018 -0.045 -0.624 -0.979
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Time              0.006   0.006     1 6562.0  0.0058   0.93939    
## Visit            62.493  62.493     1 6562.3 62.8088 2.657e-15 ***
## Age               4.898   4.898     1   24.0  4.9228   0.03621 *  
## Sex               3.682   3.682     1   24.0  3.7001   0.06635 .  
## ImageSet          2.553   2.553     1 6562.0  2.5660   0.10923    
## Time:Visit       77.782  77.782     1 6562.0 78.1749 < 2.2e-16 ***
## Age:Sex           5.193   5.193     1   24.0  5.2193   0.03147 *  
## Age:ImageSet      2.593   2.593     1 6562.0  2.6062   0.10649    
## Sex:ImageSet      4.886   4.886     1 6562.0  4.9106   0.02673 *  
## Age:Sex:ImageSet  5.279   5.279     1 6562.0  5.3056   0.02129 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug(subject_data[subject_data$ImageSet %in% c('meth', 'opioid'),], 'arousal')

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex * ImageSet)^2 + (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 20211.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.0751 -0.3333  0.0039  0.4164  7.0093 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 4.586    2.141   
##  Residual             1.206    1.098   
## Number of obs: 6600, groups:  id, 28
## 
## Fixed effects:
##                              Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)                 3.884e+00  2.526e+00  2.406e+01   1.538  0.13714   
## Time                       -6.623e-04  3.556e-04  6.565e+03  -1.863  0.06256 . 
## VisitV2                    -1.403e-01  5.441e-02  6.565e+03  -2.578  0.00995 **
## Age                        -1.766e-02  6.515e-02  2.405e+01  -0.271  0.78866   
## SexMale                    -2.643e+00  4.016e+00  2.405e+01  -0.658  0.51677   
## ImageSetopioid              1.941e-02  1.671e-01  6.565e+03   0.116  0.90751   
## Time:VisitV2                8.060e-04  5.098e-04  6.565e+03   1.581  0.11393   
## Age:SexMale                 1.414e-01  1.042e-01  2.405e+01   1.357  0.18747   
## Age:ImageSetopioid          1.039e-03  4.310e-03  6.565e+03   0.241  0.80946   
## SexMale:ImageSetopioid      3.256e-01  2.665e-01  6.565e+03   1.222  0.22190   
## Age:SexMale:ImageSetopioid -1.157e-02  6.906e-03  6.565e+03  -1.675  0.09400 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Time   VistV2 Age    SexMal ImgStp Tm:VV2 Ag:SxM Ag:ImS SxM:IS
## Time        -0.013                                                               
## VisitV2     -0.010  0.602                                                        
## Age         -0.969  0.000 -0.001                                                 
## SexMale     -0.629  0.000  0.000  0.610                                          
## ImageSetopd -0.033  0.000 -0.013  0.032  0.021                                   
## Time:VistV2  0.008 -0.698 -0.865  0.001  0.000  0.015                            
## Age:SexMale  0.606  0.000  0.000 -0.625 -0.979 -0.020 -0.001                     
## Age:ImgStpd  0.032  0.003  0.019 -0.033 -0.020 -0.970 -0.022  0.021              
## SxMl:ImgStp  0.021  0.005  0.012 -0.020 -0.033 -0.627 -0.014  0.032  0.608       
## Ag:SxMl:ImS -0.020 -0.006 -0.016  0.021  0.032  0.605  0.018 -0.033 -0.624 -0.979
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Time             1.2487  1.2487     1 6565.0  1.0354 0.308928   
## Visit            8.0167  8.0167     1 6565.2  6.6471 0.009953 **
## Age              1.1415  1.1415     1   24.0  0.9465 0.340326   
## Sex              0.4604  0.4604     1   24.0  0.3817 0.542509   
## ImageSet         2.2546  2.2546     1 6565.0  1.8694 0.171590   
## Time:Visit       3.0144  3.0144     1 6565.0  2.4994 0.113935   
## Age:Sex          2.0443  2.0443     1   24.0  1.6950 0.205298   
## Age:ImageSet     2.2771  2.2771     1 6565.0  1.8880 0.169469   
## Sex:ImageSet     1.7998  1.7998     1 6565.0  1.4923 0.221897   
## Age:Sex:ImageSet 3.3833  3.3833     1 6565.0  2.8053 0.094003 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug(subject_data[subject_data$ImageSet %in% c('meth', 'opioid'),], 'typicality')

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex * ImageSet)^2 + (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 52301.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.9362 -0.1409  0.1529  0.4851  3.6520 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 139      11.79   
##  Residual             158      12.57   
## Number of obs: 6600, groups:  id, 28
## 
## Fixed effects:
##                              Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                 7.726e+01  1.396e+01  2.427e+01   5.533 1.04e-05 ***
## Time                       -2.308e-02  4.071e-03  6.565e+03  -5.670 1.49e-08 ***
## VisitV2                    -2.251e+00  6.229e-01  6.566e+03  -3.614 0.000304 ***
## Age                         4.869e-01  3.600e-01  2.423e+01   1.352 0.188721    
## SexMale                     3.476e+00  2.219e+01  2.423e+01   0.157 0.876857    
## ImageSetopioid              6.289e-01  1.912e+00  6.565e+03   0.329 0.742270    
## Time:VisitV2               -2.672e-02  5.836e-03  6.565e+03  -4.579 4.77e-06 ***
## Age:SexMale                -1.764e-01  5.760e-01  2.423e+01  -0.306 0.761983    
## Age:ImageSetopioid          4.352e-03  4.934e-02  6.565e+03   0.088 0.929714    
## SexMale:ImageSetopioid     -2.281e+00  3.051e+00  6.565e+03  -0.748 0.454759    
## Age:SexMale:ImageSetopioid  3.339e-02  7.906e-02  6.565e+03   0.422 0.672797    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Time   VistV2 Age    SexMal ImgStp Tm:VV2 Ag:SxM Ag:ImS SxM:IS
## Time        -0.027                                                               
## VisitV2     -0.020  0.602                                                        
## Age         -0.969  0.000 -0.002                                                 
## SexMale     -0.629  0.000  0.000  0.610                                          
## ImageSetopd -0.069  0.000 -0.013  0.066  0.043                                   
## Time:VistV2  0.017 -0.698 -0.865  0.002  0.001  0.015                            
## Age:SexMale  0.606  0.000  0.001 -0.625 -0.979 -0.042 -0.001                     
## Age:ImgStpd  0.066  0.003  0.019 -0.069 -0.042 -0.970 -0.022  0.043              
## SxMl:ImgStp  0.043  0.005  0.012 -0.042 -0.069 -0.627 -0.014  0.067  0.608       
## Ag:SxMl:ImS -0.041 -0.006 -0.016  0.043  0.067  0.605  0.018 -0.069 -0.624 -0.979
## Type III Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## Time             24654.7 24654.7     1 6565.1 155.9937 < 2.2e-16 ***
## Visit             2064.0  2064.0     1 6565.8  13.0590 0.0003041 ***
## Age                320.6   320.6     1   24.0   2.0286 0.1672316    
## Sex                  1.8     1.8     1   24.0   0.0111 0.9168811    
## ImageSet            17.8    17.8     1 6565.0   0.1124 0.7374118    
## Time:Visit        3313.3  3313.3     1 6565.1  20.9635 4.767e-06 ***
## Age:Sex             12.2    12.2     1   24.0   0.0773 0.7834056    
## Age:ImageSet        44.8    44.8     1 6565.0   0.2835 0.5944163    
## Sex:ImageSet        88.3    88.3     1 6565.0   0.5588 0.4547588    
## Age:Sex:ImageSet    28.2    28.2     1 6565.0   0.1784 0.6727975    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug(subject_data[subject_data$ImageSet %in% c('meth', 'opioid'),], 'craving')

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex * ImageSet)^2 + (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 53717
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.2791 -0.1078  0.1569  0.4930  3.2842 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 130.9    11.44   
##  Residual             196.1    14.00   
## Number of obs: 6600, groups:  id, 28
## 
## Fixed effects:
##                              Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                 7.594e+01  1.357e+01  2.435e+01   5.596 8.82e-06 ***
## Time                       -2.846e-02  4.535e-03  6.565e+03  -6.276 3.69e-10 ***
## VisitV2                    -1.413e+00  6.939e-01  6.566e+03  -2.037   0.0417 *  
## Age                         4.763e-01  3.499e-01  2.430e+01   1.361   0.1859    
## SexMale                     1.086e+00  2.157e+01  2.430e+01   0.050   0.9603    
## ImageSetopioid              5.836e-01  2.130e+00  6.565e+03   0.274   0.7841    
## Time:VisitV2               -2.851e-02  6.501e-03  6.565e+03  -4.386 1.17e-05 ***
## Age:SexMale                -1.001e-01  5.597e-01  2.430e+01  -0.179   0.8596    
## Age:ImageSetopioid          1.003e-02  5.496e-02  6.565e+03   0.183   0.8552    
## SexMale:ImageSetopioid     -3.359e+00  3.399e+00  6.565e+03  -0.988   0.3231    
## Age:SexMale:ImageSetopioid  5.819e-02  8.807e-02  6.565e+03   0.661   0.5088    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Time   VistV2 Age    SexMal ImgStp Tm:VV2 Ag:SxM Ag:ImS SxM:IS
## Time        -0.031                                                               
## VisitV2     -0.023  0.602                                                        
## Age         -0.969  0.000 -0.002                                                 
## SexMale     -0.628  0.000  0.000  0.610                                          
## ImageSetopd -0.079  0.000 -0.013  0.076  0.049                                   
## Time:VistV2  0.019 -0.698 -0.865  0.003  0.001  0.015                            
## Age:SexMale  0.606  0.000  0.001 -0.625 -0.979 -0.048 -0.001                     
## Age:ImgStpd  0.076  0.003  0.019 -0.079 -0.048 -0.970 -0.022  0.049              
## SxMl:ImgStp  0.049  0.005  0.012 -0.048 -0.079 -0.627 -0.014  0.077  0.608       
## Ag:SxMl:ImS -0.047 -0.006 -0.016  0.049  0.077  0.605  0.018 -0.079 -0.624 -0.979
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## Time              33878   33878     1 6565.1 172.7388 < 2.2e-16 ***
## Visit               814     814     1 6566.0   4.1480   0.04172 *  
## Age                 501     501     1   24.0   2.5536   0.12313    
## Sex                   0       0     1   24.0   0.0008   0.97822    
## ImageSet             82      82     1 6565.0   0.4158   0.51909    
## Time:Visit         3773    3773     1 6565.1  19.2366 1.173e-05 ***
## Age:Sex               3       3     1   24.0   0.0162   0.89985    
## Age:ImageSet        155     155     1 6565.0   0.7897   0.37423    
## Sex:ImageSet        192     192     1 6565.0   0.9765   0.32309    
## Age:Sex:ImageSet     86      86     1 6565.0   0.4366   0.50879    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Neutral Cues

summarize_one_lme_neutral(subject_data[subject_data$ImageSet %in% c('control'),], 'valence')

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex) + (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 9141.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.2538 -0.1570  0.0124  0.0963  6.5360 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 1.0948   1.0463  
##  Residual             0.8862   0.9414  
## Number of obs: 3300, groups:  id, 28
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)   4.347e+00  1.237e+00  2.406e+01   3.513  0.00178 **
## Time          1.159e-03  4.728e-04  3.269e+03   2.452  0.01426 * 
## VisitV2       3.257e-02  6.516e-02  3.269e+03   0.500  0.61719   
## Age           1.340e-02  3.190e-02  2.400e+01   0.420  0.67822   
## SexMale       1.661e+00  1.967e+00  2.400e+01   0.845  0.40664   
## Time:VisitV2 -8.366e-04  6.698e-04  3.269e+03  -1.249  0.21176   
## Age:SexMale  -4.637e-02  5.104e-02  2.400e+01  -0.909  0.37264   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Time   VistV2 Age    SexMal Tm:VV2
## Time        -0.032                                   
## VisitV2     -0.028  0.612                            
## Age         -0.969  0.000  0.002                     
## SexMale     -0.628 -0.001  0.001  0.610              
## Time:VistV2  0.025 -0.706 -0.862 -0.002  0.000       
## Age:SexMale  0.606  0.000 -0.001 -0.625 -0.979  0.000
## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Time       4.3331  4.3331     1 3269.4  4.8898 0.02709 *
## Visit      0.2214  0.2214     1 3269.4  0.2499 0.61719  
## Age        0.1303  0.1303     1   24.0  0.1470 0.70475  
## Sex        0.6322  0.6322     1   24.0  0.7135 0.40664  
## Time:Visit 1.3824  1.3824     1 3269.3  1.5600 0.21176  
## Age:Sex    0.7314  0.7314     1   24.0  0.8254 0.37264  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_neutral(subject_data[subject_data$ImageSet %in% c('control'),], 'arousal')

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex) + (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 11260.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9989 -0.3580 -0.0204  0.2142  5.9661 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 1.828    1.352   
##  Residual             1.688    1.299   
## Number of obs: 3300, groups:  id, 28
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)   2.039e+00  1.600e+00  2.403e+01   1.275  0.21467   
## Time         -8.255e-05  6.525e-04  3.269e+03  -0.127  0.89933   
## VisitV2      -2.582e-02  8.993e-02  3.269e+03  -0.287  0.77405   
## Age          -1.185e-02  4.124e-02  2.396e+01  -0.287  0.77633   
## SexMale       2.685e+00  2.543e+00  2.397e+01   1.056  0.30153   
## Time:VisitV2  2.666e-03  9.244e-04  3.269e+03   2.884  0.00395 **
## Age:SexMale  -3.453e-02  6.598e-02  2.397e+01  -0.523  0.60554   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Time   VistV2 Age    SexMal Tm:VV2
## Time        -0.034                                   
## VisitV2     -0.030  0.612                            
## Age         -0.969  0.000  0.002                     
## SexMale     -0.628 -0.001  0.001  0.610              
## Time:VistV2  0.026 -0.706 -0.862 -0.002  0.000       
## Age:SexMale  0.606  0.001 -0.001 -0.625 -0.979  0.000
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Time       12.3426 12.3426     1 3269.4  7.3127 0.006882 **
## Visit       0.1391  0.1391     1 3269.4  0.0824 0.774055   
## Age         1.3146  1.3146     1   24.0  0.7789 0.386259   
## Sex         1.8819  1.8819     1   24.0  1.1150 0.301526   
## Time:Visit 14.0390 14.0390     1 3269.3  8.3178 0.003951 **
## Age:Sex     0.4623  0.4623     1   24.0  0.2739 0.605545   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_neutral(subject_data[subject_data$ImageSet %in% c('control'),], 'typicality')

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex) + (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 27913.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5429 -0.3959 -0.1737  0.0319  6.0613 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 179.3    13.39   
##  Residual             266.2    16.31   
## Number of obs: 3300, groups:  id, 28
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   3.569e+01  1.589e+01  2.409e+01   2.247 0.034098 *  
## Time         -1.387e-02  8.194e-03  3.270e+03  -1.693 0.090569 .  
## VisitV2      -4.209e+00  1.129e+00  3.270e+03  -3.727 0.000197 ***
## Age          -6.669e-01  4.094e-01  2.399e+01  -1.629 0.116328    
## SexMale      -1.732e+01  2.524e+01  2.399e+01  -0.686 0.499188    
## Time:VisitV2  5.767e-02  1.161e-02  3.270e+03   4.968 7.12e-07 ***
## Age:SexMale   5.726e-01  6.549e-01  2.399e+01   0.874 0.390600    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Time   VistV2 Age    SexMal Tm:VV2
## Time        -0.043                                   
## VisitV2     -0.038  0.612                            
## Age         -0.968  0.000  0.002                     
## SexMale     -0.628 -0.001  0.001  0.610              
## Time:VistV2  0.033 -0.706 -0.862 -0.003  0.000       
## Age:SexMale  0.605  0.001 -0.001 -0.625 -0.979  0.000
## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Time       1767.2  1767.2     1 3269.6  6.6393 0.0100189 *  
## Visit      3697.0  3697.0     1 3269.6 13.8894 0.0001972 ***
## Age         359.6   359.6     1   24.0  1.3510 0.2565287    
## Sex         125.3   125.3     1   24.0  0.4708 0.4991875    
## Time:Visit 6568.8  6568.8     1 3269.6 24.6784 7.122e-07 ***
## Age:Sex     203.5   203.5     1   24.0  0.7645 0.3905997    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_neutral(subject_data[subject_data$ImageSet %in% c('control'),], 'craving')

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + (Age * Sex) + (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 28253.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6918 -0.3422 -0.1704 -0.0075  5.7139 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 197.7    14.06   
##  Residual             295.1    17.18   
## Number of obs: 3300, groups:  id, 28
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)   2.682e+01  1.669e+01  2.409e+01   1.607  0.12102   
## Time         -9.995e-03  8.627e-03  3.270e+03  -1.159  0.24670   
## VisitV2      -2.631e+00  1.189e+00  3.270e+03  -2.212  0.02700 * 
## Age          -4.613e-01  4.299e-01  2.399e+01  -1.073  0.29397   
## SexMale      -1.198e+01  2.651e+01  2.399e+01  -0.452  0.65528   
## Time:VisitV2  3.454e-02  1.222e-02  3.270e+03   2.826  0.00475 **
## Age:SexMale   4.437e-01  6.879e-01  2.399e+01   0.645  0.52506   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Time   VistV2 Age    SexMal Tm:VV2
## Time        -0.043                                   
## VisitV2     -0.038  0.612                            
## Age         -0.968  0.000  0.002                     
## SexMale     -0.628 -0.001  0.001  0.610              
## Time:VistV2  0.033 -0.706 -0.862 -0.003  0.000       
## Age:SexMale  0.605  0.001 -0.001 -0.625 -0.979  0.000
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Time        417.58  417.58     1 3269.6  1.4151 0.234296   
## Visit      1444.44 1444.44     1 3269.6  4.8950 0.027004 * 
## Age         143.07  143.07     1   24.0  0.4848 0.492929   
## Sex          60.30   60.30     1   24.0  0.2044 0.655284   
## Time:Visit 2356.19 2356.19     1 3269.6  7.9847 0.004746 **
## Age:Sex     122.76  122.76     1   24.0  0.4160 0.525059   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Subcategories within Meth/Opioid and Injection Use

summarize_one_lme_drug_category <- function(dset, RatingType){
  this_lme <- lmer(Value ~ Time * Visit + Age * Sex + InjectionUser * Category + (1 | id), 
                data = dset[dset$RatingType == RatingType,])
  print(plot_model(this_lme, title = paste(dset$ImageSet[1], RatingType), show.values = TRUE, type = 'std'))
  print(summary(this_lme))
  #for using nlme instead of lme4--but can't get standardized betas in the forest plot
  #that_lme <- lme(fixed = Value ~ Time * Visit + Age + Sex + ImageSet, random = ~ 1 | id, 
  #              data = dset[dset$RatingType == RatingType,])
  #print(summary(that_lme))
  print(anova(this_lme))
}
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('meth'),], 'valence')
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category +      (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 9492.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.8790 -0.4014  0.0437  0.3901  6.7604 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 1.7891   1.3376  
##  Residual             0.9795   0.9897  
## Number of obs: 3297, groups:  id, 28
## 
## Fixed effects:
##                                                  Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                                     3.274e+00  1.599e+00  2.314e+01   2.048  0.05212
## Time                                            1.955e-03  4.525e-04  3.256e+03   4.320 1.61e-05
## VisitV2                                         4.048e-01  6.986e-02  3.256e+03   5.793 7.55e-09
## Age                                            -1.740e-02  4.121e-02  2.300e+01  -0.422  0.67668
## SexMale                                        -7.375e+00  2.657e+00  2.300e+01  -2.775  0.01076
## InjectionUserTRUE                               1.532e+00  6.288e-01  2.394e+01   2.436  0.02268
## Categorymeth_face_activities                    6.199e-02  1.183e-01  3.256e+03   0.524  0.60033
## Categorymeth_hand                               2.668e-01  1.366e-01  3.256e+03   1.953  0.05087
## Categorymeth_injection_hand                    -1.199e-01  1.366e-01  3.256e+03  -0.878  0.38014
## Categorymeth_instrument                         1.384e-01  1.044e-01  3.256e+03   1.325  0.18518
## Categorymeth_instrument_hand                    6.609e-03  1.183e-01  3.256e+03   0.056  0.95546
## Time:VisitV2                                   -4.202e-03  6.526e-04  3.256e+03  -6.439 1.38e-10
## Age:SexMale                                     2.179e-01  6.911e-02  2.300e+01   3.153  0.00445
## InjectionUserTRUE:Categorymeth_face_activities -6.608e-02  1.370e-01  3.256e+03  -0.482  0.62971
## InjectionUserTRUE:Categorymeth_hand            -2.483e-01  1.583e-01  3.256e+03  -1.569  0.11684
## InjectionUserTRUE:Categorymeth_injection_hand  -2.217e-01  1.583e-01  3.256e+03  -1.401  0.16129
## InjectionUserTRUE:Categorymeth_instrument      -1.966e-01  1.209e-01  3.256e+03  -1.626  0.10408
## InjectionUserTRUE:Categorymeth_instrument_hand -3.380e-02  1.370e-01  3.256e+03  -0.247  0.80517
##                                                   
## (Intercept)                                    .  
## Time                                           ***
## VisitV2                                        ***
## Age                                               
## SexMale                                        *  
## InjectionUserTRUE                              *  
## Categorymeth_face_activities                      
## Categorymeth_hand                              .  
## Categorymeth_injection_hand                       
## Categorymeth_instrument                           
## Categorymeth_instrument_hand                      
## Time:VisitV2                                   ***
## Age:SexMale                                    ** 
## InjectionUserTRUE:Categorymeth_face_activities    
## InjectionUserTRUE:Categorymeth_hand               
## InjectionUserTRUE:Categorymeth_injection_hand     
## InjectionUserTRUE:Categorymeth_instrument         
## InjectionUserTRUE:Categorymeth_instrument_hand    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE)  or
##     vcov(summary(this_lme))        if you need it

## Type III Analysis of Variance Table with Satterthwaite's method
##                        Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Time                    0.197   0.197     1 3256.2  0.2008  0.654106    
## Visit                  32.875  32.875     1 3256.4 33.5641 7.552e-09 ***
## Age                     7.529   7.529     1   23.0  7.6871  0.010831 *  
## Sex                     7.545   7.545     1   23.0  7.7028  0.010760 *  
## InjectionUser           4.978   4.978     1   23.0  5.0823  0.033986 *  
## Category               19.930   3.986     5 3256.0  4.0695  0.001101 ** 
## Time:Visit             40.613  40.613     1 3256.1 41.4637 1.377e-10 ***
## Age:Sex                 9.736   9.736     1   23.0  9.9405  0.004452 ** 
## InjectionUser:Category  5.352   1.070     5 3256.0  1.0927  0.362171    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('meth'),], 'arousal')
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category +      (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 10149.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.3380 -0.3677  0.0272  0.3884  6.8963 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 4.323    2.079   
##  Residual             1.188    1.090   
## Number of obs: 3300, groups:  id, 28
## 
## Fixed effects:
##                                                  Estimate Std. Error         df t value Pr(>|t|)  
## (Intercept)                                     3.220e+00  2.481e+00  2.307e+01   1.298   0.2072  
## Time                                           -2.886e-04  4.981e-04  3.259e+03  -0.579   0.5623  
## VisitV2                                        -7.477e-02  7.691e-02  3.259e+03  -0.972   0.3310  
## Age                                            -3.242e-02  6.398e-02  2.300e+01  -0.507   0.6172  
## SexMale                                        -4.722e+00  4.126e+00  2.300e+01  -1.145   0.2642  
## InjectionUserTRUE                               1.637e+00  9.716e-01  2.347e+01   1.685   0.1052  
## Categorymeth_face_activities                    1.790e-01  1.303e-01  3.259e+03   1.374   0.1694  
## Categorymeth_hand                               1.521e-01  1.504e-01  3.259e+03   1.011   0.3121  
## Categorymeth_injection_hand                    -1.696e-02  1.504e-01  3.259e+03  -0.113   0.9102  
## Categorymeth_instrument                        -3.046e-02  1.149e-01  3.259e+03  -0.265   0.7910  
## Categorymeth_instrument_hand                    7.707e-02  1.303e-01  3.259e+03   0.592   0.5542  
## Time:VisitV2                                    5.338e-04  7.183e-04  3.259e+03   0.743   0.4574  
## Age:SexMale                                     1.965e-01  1.073e-01  2.300e+01   1.831   0.0801 .
## InjectionUserTRUE:Categorymeth_face_activities -3.047e-01  1.509e-01  3.259e+03  -2.020   0.0435 *
## InjectionUserTRUE:Categorymeth_hand            -2.955e-01  1.743e-01  3.259e+03  -1.696   0.0900 .
## InjectionUserTRUE:Categorymeth_injection_hand   1.375e-01  1.743e-01  3.259e+03   0.789   0.4303  
## InjectionUserTRUE:Categorymeth_instrument      -1.439e-01  1.331e-01  3.259e+03  -1.081   0.2797  
## InjectionUserTRUE:Categorymeth_instrument_hand -2.292e-01  1.509e-01  3.259e+03  -1.519   0.1289  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE)  or
##     vcov(summary(this_lme))        if you need it

## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Time                    0.0043  0.0043     1 3259.1  0.0037 0.95182  
## Visit                   1.1225  1.1225     1 3259.2  0.9451 0.33103  
## Age                     1.9568  1.9568     1   23.0  1.6476 0.21206  
## Sex                     1.5558  1.5558     1   23.0  1.3100 0.26417  
## InjectionUser           2.8509  2.8509     1   23.0  2.4005 0.13494  
## Category                7.6126  1.5225     5 3259.0  1.2820 0.26865  
## Time:Visit              0.6560  0.6560     1 3259.1  0.5524 0.45741  
## Age:Sex                 3.9806  3.9806     1   23.0  3.3517 0.08012 .
## InjectionUser:Category 12.0005  2.4001     5 3259.0  2.0209 0.07263 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('meth'),], 'typicality')
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category +      (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 26292.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.8094 -0.1627  0.1535  0.4861  2.1803 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 152.2    12.33   
##  Residual             164.0    12.81   
## Number of obs: 3300, groups:  id, 28
## 
## Fixed effects:
##                                                  Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                                     7.795e+01  1.480e+01  2.328e+01   5.267 2.33e-05
## Time                                           -2.021e-02  5.854e-03  3.259e+03  -3.453 0.000562
## VisitV2                                        -2.037e+00  9.038e-01  3.260e+03  -2.254 0.024275
## Age                                             5.097e-01  3.808e-01  2.300e+01   1.338 0.193879
## SexMale                                         6.700e+00  2.456e+01  2.300e+01   0.273 0.787409
## InjectionUserTRUE                              -8.166e-01  5.867e+00  2.485e+01  -0.139 0.890414
## Categorymeth_face_activities                    1.212e+00  1.531e+00  3.259e+03   0.792 0.428425
## Categorymeth_hand                               1.153e+00  1.768e+00  3.259e+03   0.652 0.514150
## Categorymeth_injection_hand                     1.782e+00  1.768e+00  3.259e+03   1.008 0.313462
## Categorymeth_instrument                        -1.790e+00  1.351e+00  3.259e+03  -1.325 0.185166
## Categorymeth_instrument_hand                    7.153e-02  1.531e+00  3.259e+03   0.047 0.962741
## Time:VisitV2                                   -3.023e-02  8.441e-03  3.259e+03  -3.581 0.000348
## Age:SexMale                                    -2.617e-01  6.387e-01  2.300e+01  -0.410 0.685802
## InjectionUserTRUE:Categorymeth_face_activities -1.685e+00  1.773e+00  3.259e+03  -0.950 0.342153
## InjectionUserTRUE:Categorymeth_hand            -2.984e+00  2.048e+00  3.259e+03  -1.457 0.145144
## InjectionUserTRUE:Categorymeth_injection_hand  -3.030e+00  2.048e+00  3.259e+03  -1.479 0.139136
## InjectionUserTRUE:Categorymeth_instrument      -4.023e-01  1.564e+00  3.259e+03  -0.257 0.796991
## InjectionUserTRUE:Categorymeth_instrument_hand -2.953e+00  1.773e+00  3.259e+03  -1.665 0.095947
##                                                   
## (Intercept)                                    ***
## Time                                           ***
## VisitV2                                        *  
## Age                                               
## SexMale                                           
## InjectionUserTRUE                                 
## Categorymeth_face_activities                      
## Categorymeth_hand                                 
## Categorymeth_injection_hand                       
## Categorymeth_instrument                           
## Categorymeth_instrument_hand                      
## Time:VisitV2                                   ***
## Age:SexMale                                       
## InjectionUserTRUE:Categorymeth_face_activities    
## InjectionUserTRUE:Categorymeth_hand               
## InjectionUserTRUE:Categorymeth_injection_hand     
## InjectionUserTRUE:Categorymeth_instrument         
## InjectionUserTRUE:Categorymeth_instrument_hand .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE)  or
##     vcov(summary(this_lme))        if you need it

## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Time                   11466.4 11466.4     1 3259.3 69.9058 < 2.2e-16 ***
## Visit                    833.2   833.2     1 3259.8  5.0796 0.0242748 *  
## Age                      252.8   252.8     1   23.0  1.5409 0.2269949    
## Sex                       12.2    12.2     1   23.0  0.0744 0.7874094    
## InjectionUser             35.0    35.0     1   23.1  0.2133 0.6485393    
## Category                2389.7   477.9     5 3259.0  2.9138 0.0125026 *  
## Time:Visit              2103.2  2103.2     1 3259.2 12.8223 0.0003475 ***
## Age:Sex                   27.5    27.5     1   23.0  0.1679 0.6858018    
## InjectionUser:Category   968.8   193.8     5 3259.0  1.1813 0.3156999    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('meth'),], 'craving')
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category +      (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 26970.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.9852 -0.1233  0.1477  0.5068  3.0025 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 140.9    11.87   
##  Residual             202.1    14.22   
## Number of obs: 3300, groups:  id, 28
## 
## Fixed effects:
##                                                  Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                                     7.650e+01  1.428e+01  2.337e+01   5.358 1.83e-05
## Time                                           -2.095e-02  6.498e-03  3.259e+03  -3.225  0.00127
## VisitV2                                         1.267e-01  1.003e+00  3.260e+03   0.126  0.89953
## Age                                             4.821e-01  3.670e-01  2.300e+01   1.313  0.20198
## SexMale                                         2.004e+00  2.367e+01  2.300e+01   0.085  0.93326
## InjectionUserTRUE                               2.737e-01  5.689e+00  2.548e+01   0.048  0.96200
## Categorymeth_face_activities                    2.876e-01  1.699e+00  3.259e+03   0.169  0.86561
## Categorymeth_hand                               2.905e-02  1.962e+00  3.259e+03   0.015  0.98819
## Categorymeth_injection_hand                     2.994e-01  1.962e+00  3.259e+03   0.153  0.87874
## Categorymeth_instrument                        -3.114e+00  1.499e+00  3.259e+03  -2.077  0.03784
## Categorymeth_instrument_hand                   -1.367e+00  1.700e+00  3.259e+03  -0.805  0.42115
## Time:VisitV2                                   -4.183e-02  9.370e-03  3.259e+03  -4.465 8.29e-06
## Age:SexMale                                    -1.245e-01  6.156e-01  2.300e+01  -0.202  0.84147
## InjectionUserTRUE:Categorymeth_face_activities -2.100e+00  1.968e+00  3.259e+03  -1.067  0.28604
## InjectionUserTRUE:Categorymeth_hand            -2.694e+00  2.273e+00  3.259e+03  -1.185  0.23604
## InjectionUserTRUE:Categorymeth_injection_hand  -1.215e+00  2.273e+00  3.259e+03  -0.534  0.59318
## InjectionUserTRUE:Categorymeth_instrument       2.638e-01  1.736e+00  3.259e+03   0.152  0.87923
## InjectionUserTRUE:Categorymeth_instrument_hand -2.108e+00  1.968e+00  3.259e+03  -1.071  0.28414
##                                                   
## (Intercept)                                    ***
## Time                                           ** 
## VisitV2                                           
## Age                                               
## SexMale                                           
## InjectionUserTRUE                                 
## Categorymeth_face_activities                      
## Categorymeth_hand                                 
## Categorymeth_injection_hand                       
## Categorymeth_instrument                        *  
## Categorymeth_instrument_hand                      
## Time:VisitV2                                   ***
## Age:SexMale                                       
## InjectionUserTRUE:Categorymeth_face_activities    
## InjectionUserTRUE:Categorymeth_hand               
## InjectionUserTRUE:Categorymeth_injection_hand     
## InjectionUserTRUE:Categorymeth_instrument         
## InjectionUserTRUE:Categorymeth_instrument_hand    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE)  or
##     vcov(summary(this_lme))        if you need it

## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Time                   16109.0 16109.0     1 3259.4 79.7049 < 2.2e-16 ***
## Visit                      3.2     3.2     1 3260.1  0.0159  0.899528    
## Age                      411.8   411.8     1   23.0  2.0375  0.166898    
## Sex                        1.4     1.4     1   23.0  0.0072  0.933264    
## InjectionUser              7.0     7.0     1   23.1  0.0348  0.853668    
## Category                3508.8   701.8     5 3259.0  3.4722  0.003926 ** 
## Time:Visit              4028.7  4028.7     1 3259.3 19.9333 8.291e-06 ***
## Age:Sex                    8.3     8.3     1   23.0  0.0409  0.841467    
## InjectionUser:Category   840.6   168.1     5 3259.0  0.8318  0.526839    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('opioid'),], 'valence')
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category +      (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 9405.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0560 -0.3698  0.0620  0.4468  4.4949 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 1.8191   1.3487  
##  Residual             0.9511   0.9753  
## Number of obs: 3300, groups:  id, 28
## 
## Fixed effects:
##                                                         Estimate Std. Error         df t value
## (Intercept)                                            3.471e+00  1.612e+00  2.313e+01   2.153
## Time                                                   2.090e-03  4.501e-04  3.259e+03   4.642
## VisitV2                                                3.703e-01  6.814e-02  3.259e+03   5.435
## Age                                                   -1.310e-02  4.154e-02  2.300e+01  -0.315
## SexMale                                               -6.565e+00  2.679e+00  2.300e+01  -2.451
## InjectionUserTRUE                                      1.321e+00  6.337e-01  2.390e+01   2.084
## Categoryopioid_face_activities                        -2.766e-01  1.166e-01  3.259e+03  -2.373
## Categoryopioid_hand                                    2.285e-01  1.346e-01  3.259e+03   1.697
## Categoryopioid_injection_hand                         -3.905e-01  1.346e-01  3.259e+03  -2.901
## Categoryopioid_injection_instrument                   -1.485e-01  1.028e-01  3.259e+03  -1.444
## Categoryopioid_instrument_hand                        -1.613e-01  1.166e-01  3.259e+03  -1.384
## Time:VisitV2                                          -3.894e-03  6.415e-04  3.259e+03  -6.070
## Age:SexMale                                            1.960e-01  6.968e-02  2.300e+01   2.814
## InjectionUserTRUE:Categoryopioid_face_activities      -1.016e-01  1.350e-01  3.259e+03  -0.753
## InjectionUserTRUE:Categoryopioid_hand                 -2.469e-01  1.559e-01  3.259e+03  -1.584
## InjectionUserTRUE:Categoryopioid_injection_hand       -1.306e-01  1.559e-01  3.259e+03  -0.838
## InjectionUserTRUE:Categoryopioid_injection_instrument -1.424e-01  1.191e-01  3.259e+03  -1.196
## InjectionUserTRUE:Categoryopioid_instrument_hand      -1.607e-01  1.351e-01  3.259e+03  -1.190
##                                                       Pr(>|t|)    
## (Intercept)                                            0.04200 *  
## Time                                                  3.58e-06 ***
## VisitV2                                               5.89e-08 ***
## Age                                                    0.75535    
## SexMale                                                0.02228 *  
## InjectionUserTRUE                                      0.04803 *  
## Categoryopioid_face_activities                         0.01771 *  
## Categoryopioid_hand                                    0.08975 .  
## Categoryopioid_injection_hand                          0.00374 ** 
## Categoryopioid_injection_instrument                    0.14885    
## Categoryopioid_instrument_hand                         0.16660    
## Time:VisitV2                                          1.43e-09 ***
## Age:SexMale                                            0.00986 ** 
## InjectionUserTRUE:Categoryopioid_face_activities       0.45173    
## InjectionUserTRUE:Categoryopioid_hand                  0.11336    
## InjectionUserTRUE:Categoryopioid_injection_hand        0.40224    
## InjectionUserTRUE:Categoryopioid_injection_instrument  0.23178    
## InjectionUserTRUE:Categoryopioid_instrument_hand       0.23413    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE)  or
##     vcov(summary(this_lme))        if you need it

## Type III Analysis of Variance Table with Satterthwaite's method
##                        Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Time                    0.188   0.188     1 3259.1  0.1982  0.656224    
## Visit                  28.093  28.093     1 3259.4 29.5366 5.891e-08 ***
## Age                     6.189   6.189     1   23.0  6.5076  0.017859 *  
## Sex                     5.712   5.712     1   23.0  6.0059  0.022281 *  
## InjectionUser           3.418   3.418     1   23.0  3.5939  0.070611 .  
## Category               64.145  12.829     5 3259.0 13.4883 4.779e-13 ***
## Time:Visit             35.042  35.042     1 3259.2 36.8430 1.428e-09 ***
## Age:Sex                 7.529   7.529     1   23.0  7.9159  0.009859 ** 
## InjectionUser:Category  2.809   0.562     5 3259.0  0.5907  0.707114    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('opioid'),], 'arousal')
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category +      (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 10093.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.4731 -0.2912  0.0237  0.4270  3.5615 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 4.272    2.067   
##  Residual             1.167    1.080   
## Number of obs: 3300, groups:  id, 28
## 
## Fixed effects:
##                                                         Estimate Std. Error         df t value
## (Intercept)                                            3.271e+00  2.466e+00  2.307e+01   1.326
## Time                                                  -1.274e-03  4.987e-04  3.259e+03  -2.555
## VisitV2                                               -2.422e-01  7.550e-02  3.259e+03  -3.208
## Age                                                   -3.333e-02  6.360e-02  2.300e+01  -0.524
## SexMale                                               -4.659e+00  4.101e+00  2.300e+01  -1.136
## InjectionUserTRUE                                      1.557e+00  9.658e-01  2.347e+01   1.612
## Categoryopioid_face_activities                         1.214e-02  1.292e-01  3.259e+03   0.094
## Categoryopioid_hand                                    1.236e-01  1.491e-01  3.259e+03   0.829
## Categoryopioid_injection_hand                          3.257e-02  1.491e-01  3.259e+03   0.218
## Categoryopioid_injection_instrument                    7.458e-02  1.139e-01  3.259e+03   0.655
## Categoryopioid_instrument_hand                         2.097e-01  1.292e-01  3.259e+03   1.623
## Time:VisitV2                                           1.486e-03  7.107e-04  3.259e+03   2.092
## Age:SexMale                                            1.919e-01  1.067e-01  2.300e+01   1.799
## InjectionUserTRUE:Categoryopioid_face_activities       2.936e-01  1.496e-01  3.259e+03   1.963
## InjectionUserTRUE:Categoryopioid_hand                 -1.746e-01  1.727e-01  3.259e+03  -1.011
## InjectionUserTRUE:Categoryopioid_injection_hand        4.901e-01  1.727e-01  3.259e+03   2.838
## InjectionUserTRUE:Categoryopioid_injection_instrument  1.494e-01  1.319e-01  3.259e+03   1.132
## InjectionUserTRUE:Categoryopioid_instrument_hand      -1.336e-03  1.496e-01  3.259e+03  -0.009
##                                                       Pr(>|t|)   
## (Intercept)                                            0.19775   
## Time                                                   0.01066 * 
## VisitV2                                                0.00135 **
## Age                                                    0.60527   
## SexMale                                                0.26765   
## InjectionUserTRUE                                      0.12039   
## Categoryopioid_face_activities                         0.92511   
## Categoryopioid_hand                                    0.40723   
## Categoryopioid_injection_hand                          0.82711   
## Categoryopioid_injection_instrument                    0.51274   
## Categoryopioid_instrument_hand                         0.10474   
## Time:VisitV2                                           0.03656 * 
## Age:SexMale                                            0.08515 . 
## InjectionUserTRUE:Categoryopioid_face_activities       0.04977 * 
## InjectionUserTRUE:Categoryopioid_hand                  0.31206   
## InjectionUserTRUE:Categoryopioid_injection_hand        0.00457 **
## InjectionUserTRUE:Categoryopioid_injection_instrument  0.25759   
## InjectionUserTRUE:Categoryopioid_instrument_hand       0.99288   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE)  or
##     vcov(summary(this_lme))        if you need it

## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Time                    2.6064  2.6064     1 3259.1  2.2326 0.135222   
## Visit                  12.0177 12.0177     1 3259.2 10.2943 0.001347 **
## Age                     1.7628  1.7628     1   23.0  1.5100 0.231558   
## Sex                     1.5066  1.5066     1   23.0  1.2905 0.267647   
## InjectionUser           3.5789  3.5789     1   23.0  3.0657 0.093281 . 
## Category               17.8034  3.5607     5 3259.0  3.0501 0.009456 **
## Time:Visit              5.1069  5.1069     1 3259.1  4.3745 0.036557 * 
## Age:Sex                 3.7784  3.7784     1   23.0  3.2365 0.085149 . 
## InjectionUser:Category 20.3376  4.0675     5 3259.0  3.4842 0.003828 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('opioid'),], 'typicality')
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category +      (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 25997.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.8348 -0.1406  0.1551  0.4661  3.8109 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 136.3    11.67   
##  Residual             150.0    12.25   
## Number of obs: 3300, groups:  id, 28
## 
## Fixed effects:
##                                                         Estimate Std. Error         df t value
## (Intercept)                                            7.898e+01  1.401e+01  2.328e+01   5.638
## Time                                                  -2.598e-02  5.652e-03  3.259e+03  -4.596
## VisitV2                                               -2.593e+00  8.556e-01  3.260e+03  -3.030
## Age                                                    5.158e-01  3.604e-01  2.300e+01   1.431
## SexMale                                                4.633e+00  2.324e+01  2.300e+01   0.199
## InjectionUserTRUE                                     -4.029e+00  5.555e+00  2.489e+01  -0.725
## Categoryopioid_face_activities                         1.616e+00  1.464e+00  3.259e+03   1.104
## Categoryopioid_hand                                   -1.970e+00  1.690e+00  3.259e+03  -1.165
## Categoryopioid_injection_hand                          1.826e+00  1.690e+00  3.259e+03   1.080
## Categoryopioid_injection_instrument                   -1.261e+00  1.291e+00  3.259e+03  -0.977
## Categoryopioid_instrument_hand                         1.095e+00  1.464e+00  3.259e+03   0.748
## Time:VisitV2                                          -2.181e-02  8.055e-03  3.259e+03  -2.707
## Age:SexMale                                           -2.342e-01  6.045e-01  2.300e+01  -0.387
## InjectionUserTRUE:Categoryopioid_face_activities       7.640e-01  1.695e+00  3.259e+03   0.451
## InjectionUserTRUE:Categoryopioid_hand                  8.393e-01  1.958e+00  3.259e+03   0.429
## InjectionUserTRUE:Categoryopioid_injection_hand        1.283e+00  1.958e+00  3.259e+03   0.655
## InjectionUserTRUE:Categoryopioid_injection_instrument  3.275e+00  1.495e+00  3.259e+03   2.190
## InjectionUserTRUE:Categoryopioid_instrument_hand       1.285e+00  1.696e+00  3.259e+03   0.758
##                                                       Pr(>|t|)    
## (Intercept)                                           9.31e-06 ***
## Time                                                  4.48e-06 ***
## VisitV2                                                0.00246 ** 
## Age                                                    0.16583    
## SexMale                                                0.84374    
## InjectionUserTRUE                                      0.47502    
## Categoryopioid_face_activities                         0.26977    
## Categoryopioid_hand                                    0.24393    
## Categoryopioid_injection_hand                          0.28015    
## Categoryopioid_injection_instrument                    0.32881    
## Categoryopioid_instrument_hand                         0.45476    
## Time:VisitV2                                           0.00682 ** 
## Age:SexMale                                            0.70205    
## InjectionUserTRUE:Categoryopioid_face_activities       0.65228    
## InjectionUserTRUE:Categoryopioid_hand                  0.66812    
## InjectionUserTRUE:Categoryopioid_injection_hand        0.51222    
## InjectionUserTRUE:Categoryopioid_injection_instrument  0.02857 *  
## InjectionUserTRUE:Categoryopioid_instrument_hand       0.44862    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE)  or
##     vcov(summary(this_lme))        if you need it

## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Time                   12568.7 12568.7     1 3259.2 83.8112 < 2.2e-16 ***
## Visit                   1377.2  1377.2     1 3259.9  9.1832  0.002462 ** 
## Age                      285.9   285.9     1   23.0  1.9061  0.180666    
## Sex                        6.0     6.0     1   23.0  0.0397  0.843737    
## InjectionUser             39.3    39.3     1   23.1  0.2617  0.613788    
## Category                3445.4   689.1     5 3259.0  4.5949  0.000352 ***
## Time:Visit              1099.2  1099.2     1 3259.3  7.3300  0.006817 ** 
## Age:Sex                   22.5    22.5     1   23.0  0.1501  0.702047    
## InjectionUser:Category   917.3   183.5     5 3259.0  1.2233  0.295304    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_drug_category(subject_data[subject_data$ImageSet %in% c('opioid'),], 'craving')
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time * Visit + Age * Sex + InjectionUser * Category +      (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 26704
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.5357 -0.1324  0.1776  0.4725  3.4332 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 133.1    11.54   
##  Residual             186.3    13.65   
## Number of obs: 3300, groups:  id, 28
## 
## Fixed effects:
##                                                         Estimate Std. Error         df t value
## (Intercept)                                            7.700e+01  1.387e+01  2.336e+01   5.551
## Time                                                  -3.645e-02  6.300e-03  3.259e+03  -5.786
## VisitV2                                               -3.165e+00  9.537e-01  3.260e+03  -3.319
## Age                                                    4.851e-01  3.566e-01  2.300e+01   1.360
## SexMale                                               -2.330e+00  2.300e+01  2.300e+01  -0.101
## InjectionUserTRUE                                     -1.829e+00  5.525e+00  2.542e+01  -0.331
## Categoryopioid_face_activities                         2.885e+00  1.632e+00  3.259e+03   1.768
## Categoryopioid_hand                                   -1.747e+00  1.884e+00  3.259e+03  -0.927
## Categoryopioid_injection_hand                          1.484e+00  1.884e+00  3.259e+03   0.788
## Categoryopioid_injection_instrument                   -1.286e+00  1.439e+00  3.259e+03  -0.894
## Categoryopioid_instrument_hand                         1.925e+00  1.632e+00  3.259e+03   1.180
## Time:VisitV2                                          -1.274e-02  8.978e-03  3.259e+03  -1.419
## Age:SexMale                                           -3.938e-02  5.981e-01  2.300e+01  -0.066
## InjectionUserTRUE:Categoryopioid_face_activities       3.055e-01  1.890e+00  3.259e+03   0.162
## InjectionUserTRUE:Categoryopioid_hand                  1.123e+00  2.182e+00  3.259e+03   0.515
## InjectionUserTRUE:Categoryopioid_injection_hand        3.472e+00  2.182e+00  3.259e+03   1.591
## InjectionUserTRUE:Categoryopioid_injection_instrument  3.516e+00  1.666e+00  3.259e+03   2.110
## InjectionUserTRUE:Categoryopioid_instrument_hand       1.891e+00  1.890e+00  3.259e+03   1.000
##                                                       Pr(>|t|)    
## (Intercept)                                           1.14e-05 ***
## Time                                                  7.87e-09 ***
## VisitV2                                               0.000914 ***
## Age                                                   0.186878    
## SexMale                                               0.920165    
## InjectionUserTRUE                                     0.743330    
## Categoryopioid_face_activities                        0.077145 .  
## Categoryopioid_hand                                   0.353876    
## Categoryopioid_injection_hand                         0.430899    
## Categoryopioid_injection_instrument                   0.371506    
## Categoryopioid_instrument_hand                        0.238175    
## Time:VisitV2                                          0.156081    
## Age:SexMale                                           0.948075    
## InjectionUserTRUE:Categoryopioid_face_activities      0.871566    
## InjectionUserTRUE:Categoryopioid_hand                 0.606749    
## InjectionUserTRUE:Categoryopioid_injection_hand       0.111678    
## InjectionUserTRUE:Categoryopioid_injection_instrument 0.034932 *  
## InjectionUserTRUE:Categoryopioid_instrument_hand      0.317242    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(summary(this_lme), correlation=TRUE)  or
##     vcov(summary(this_lme))        if you need it

## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Time                   16945.9 16945.9     1 3259.3 90.9614 < 2.2e-16 ***
## Visit                   2051.9  2051.9     1 3260.1 11.0143 0.0009141 ***
## Age                      494.3   494.3     1   23.0  2.6530 0.1169728    
## Sex                        1.9     1.9     1   23.0  0.0103 0.9201647    
## InjectionUser              0.1     0.1     1   23.1  0.0004 0.9837537    
## Category                5978.7  1195.7     5 3259.0  6.4185 6.078e-06 ***
## Time:Visit               375.0   375.0     1 3259.4  2.0127 0.1560811    
## Age:Sex                    0.8     0.8     1   23.0  0.0043 0.9480752    
## InjectionUser:Category  1324.8   265.0     5 3259.0  1.4223 0.2128160    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Image Characteristics–Hue, Saturation, Valence

summarize_one_lm_imagechar <- function(dset, characteristic){
  this_lm <- lm(as.formula(paste(characteristic, '~ ImageSet')), 
                data = dset)
  print(plot_model(this_lm, title = characteristic, show.values = TRUE, type = 'est'))
  print(summary(this_lm))
  print(anova(this_lm))
}
summarize_one_lm_imagechar(merged_summary_values, 'hue_mean')

## 
## Call:
## lm(formula = as.formula(paste(characteristic, "~ ImageSet")), 
##     data = dset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27365 -0.14219 -0.04265  0.10681  0.59849 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     0.19051    0.01493  12.757  < 2e-16 ***
## ImageSetmeth    0.08313    0.02112   3.936 9.85e-05 ***
## ImageSetopioid  0.06968    0.02112   3.299  0.00106 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1683 on 378 degrees of freedom
## Multiple R-squared:  0.04511,    Adjusted R-squared:  0.04005 
## F-statistic: 8.928 on 2 and 378 DF,  p-value: 0.0001628
## 
## Analysis of Variance Table
## 
## Response: hue_mean
##            Df  Sum Sq  Mean Sq F value    Pr(>F)    
## ImageSet    2  0.5058 0.252875  8.9277 0.0001628 ***
## Residuals 378 10.7068 0.028325                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lm_imagechar(merged_summary_values, 'saturation_mean')

## 
## Call:
## lm(formula = as.formula(paste(characteristic, "~ ImageSet")), 
##     data = dset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30227 -0.12424 -0.01327  0.10830  0.44876 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     0.18970    0.01446  13.117  < 2e-16 ***
## ImageSetmeth    0.11257    0.02045   5.504 6.86e-08 ***
## ImageSetopioid  0.10954    0.02045   5.356 1.48e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.163 on 378 degrees of freedom
## Multiple R-squared:  0.09425,    Adjusted R-squared:  0.08946 
## F-statistic: 19.67 on 2 and 378 DF,  p-value: 7.493e-09
## 
## Analysis of Variance Table
## 
## Response: saturation_mean
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## ImageSet    2  1.0448 0.52240  19.667 7.493e-09 ***
## Residuals 378 10.0406 0.02656                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lm_imagechar(merged_summary_values, 'value_mean')

## 
## Call:
## lm(formula = as.formula(paste(characteristic, "~ ImageSet")), 
##     data = dset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42727 -0.15408  0.00472  0.10972  0.59092 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     0.45227    0.01884  24.004   <2e-16 ***
## ImageSetmeth   -0.05318    0.02665  -1.996   0.0467 *  
## ImageSetopioid -0.02398    0.02665  -0.900   0.3687    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2123 on 378 degrees of freedom
## Multiple R-squared:  0.01046,    Adjusted R-squared:  0.005226 
## F-statistic: 1.998 on 2 and 378 DF,  p-value: 0.137
## 
## Analysis of Variance Table
## 
## Response: value_mean
##            Df  Sum Sq  Mean Sq F value Pr(>F)
## ImageSet    2  0.1802 0.090085  1.9982  0.137
## Residuals 378 17.0419 0.045084
#ddq_data <- read.csv('DDQ_Scores.csv')
ddq_data <- read.csv('DCR_DDQScores-anonymized.csv')

summarize_one_lme_ddq <- function(dset, subscore){
  this_lme <- lmer(Value ~ Time + (1 | id), 
                data = dset[dset$Subscore == subscore,])
  print(plot_model(this_lme, title = subscore, show.values = TRUE, type = 'std'))
  print(summary(this_lme))
  print(anova(this_lme))
}

summarize_one_lme_ddq(ddq_data, 'Desire')

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time + (1 | id)
##    Data: dset[dset$Subscore == subscore, ]
## 
## REML criterion at convergence: 215.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.94412 -0.25263 -0.06047  0.21956  3.10401 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.6939   0.8330  
##  Residual             0.1921   0.4383  
## Number of obs: 112, groups:  id, 28
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  1.70918    0.17789 38.02770   9.608 1.02e-11 ***
## TimeT2       0.03571    0.11715 81.00000   0.305    0.761    
## TimeT3      -0.07143    0.11715 81.00000  -0.610    0.544    
## TimeT4      -0.08673    0.11715 81.00000  -0.740    0.461    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) TimeT2 TimeT3
## TimeT2 -0.329              
## TimeT3 -0.329  0.500       
## TimeT4 -0.329  0.500  0.500
## Type III Analysis of Variance Table with Satterthwaite's method
##       Sum Sq  Mean Sq NumDF DenDF F value Pr(>F)
## Time 0.28426 0.094752     3    81  0.4932  0.688
summarize_one_lme_ddq(ddq_data, 'Neg')

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time + (1 | id)
##    Data: dset[dset$Subscore == subscore, ]
## 
## REML criterion at convergence: 338.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3745 -0.3277 -0.0669  0.2760  2.3180 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 2.5480   1.5963  
##  Residual             0.5674   0.7533  
## Number of obs: 112, groups:  id, 28
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   2.3929     0.3336 35.9197   7.174 1.97e-08 ***
## TimeT2       -0.1339     0.2013 81.0000  -0.665   0.5078    
## TimeT3       -0.1964     0.2013 81.0000  -0.976   0.3321    
## TimeT4       -0.4018     0.2013 81.0000  -1.996   0.0493 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) TimeT2 TimeT3
## TimeT2 -0.302              
## TimeT3 -0.302  0.500       
## TimeT4 -0.302  0.500  0.500
## Type III Analysis of Variance Table with Satterthwaite's method
##      Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 2.3504 0.78348     3    81  1.3807 0.2546
summarize_one_lme_ddq(ddq_data, 'Control')

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ Time + (1 | id)
##    Data: dset[dset$Subscore == subscore, ]
## 
## REML criterion at convergence: 332.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.06379 -0.36536 -0.07568  0.27796  2.94957 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 2.3385   1.5292  
##  Residual             0.5472   0.7397  
## Number of obs: 112, groups:  id, 28
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  2.14286    0.32103 36.36237   6.675 8.38e-08 ***
## TimeT2       0.16071    0.19770 81.00000   0.813    0.419    
## TimeT3      -0.05357    0.19770 81.00000  -0.271    0.787    
## TimeT4      -0.30357    0.19770 81.00000  -1.535    0.129    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) TimeT2 TimeT3
## TimeT2 -0.308              
## TimeT3 -0.308  0.500       
## TimeT4 -0.308  0.500  0.500
## Type III Analysis of Variance Table with Satterthwaite's method
##      Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 3.1138  1.0379     3    81  1.8968 0.1367

new LMEs for 7/2019

need: all images craving|arousal|valence|absvalence ~ drug (M/O/N) + age + gender

library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:sm':
## 
##     muscle
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
## The following object is masked from 'package:sm':
## 
##     geyser
summarize_one_lme_allimages <- function(dset, RatingType){
  this_lme <- lmer(Value ~  ImageSet + Age + Sex + (1 | id), 
                data = dset[dset$RatingType == RatingType,])
  print(paste("LME OF ALL IMAGES FOR:", RatingType))
  print(plot_model(this_lme, title = paste(RatingType), show.values = TRUE, type = 'std', order.terms = c(3,4,1,2)))
  print(summary(this_lme))
  print(anova(this_lme))
  
  this_lme.ht <- glht(this_lme, linfct = mcp(ImageSet = "Tukey"))
  summary(this_lme.ht, test=adjusted("none"))
  summary(this_lme.ht)
}
summarize_one_lme_allimages(subject_data, 'craving')
## [1] "LME OF ALL IMAGES FOR: craving"

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ ImageSet + Age + Sex + (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 86227.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0716 -0.6391  0.2088  0.4583  5.3166 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept)  32.66    5.715  
##  Residual             351.75   18.755  
## Number of obs: 9900, groups:  id, 28
## 
## Fixed effects:
##                 Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)       3.9670     5.4414   25.1223   0.729    0.473    
## ImageSetmeth     76.1436     0.4617 9870.0491 164.914   <2e-16 ***
## ImageSetopioid   76.4842     0.4617 9870.0491 165.652   <2e-16 ***
## Age               0.2087     0.1376   25.0025   1.517    0.142    
## SexMale          -0.4920     2.2159   25.0336  -0.222    0.826    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ImgStm ImgStp Age   
## ImageSetmth -0.042                     
## ImageSetopd -0.042  0.500              
## Age         -0.950  0.000  0.000       
## SexMale     -0.219  0.000  0.000 -0.015
## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF DenDF    F value Pr(>F)    
## ImageSet 12812590 6406295     2  9870 18212.5694 <2e-16 ***
## Age           809     809     1    25     2.3011 0.1418    
## Sex            17      17     1    25     0.0493 0.8261    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = Value ~ ImageSet + Age + Sex + (1 | id), data = dset[dset$RatingType == 
##     RatingType, ])
## 
## Linear Hypotheses:
##                       Estimate Std. Error z value Pr(>|z|)    
## meth - control == 0    76.1436     0.4617 164.914   <1e-04 ***
## opioid - control == 0  76.4842     0.4617 165.652   <1e-04 ***
## opioid - meth == 0      0.3406     0.4617   0.738    0.741    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
summarize_one_lme_allimages(subject_data, 'arousal')
## [1] "LME OF ALL IMAGES FOR: arousal"

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ ImageSet + Age + Sex + (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 36788.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9768 -0.5700  0.0353  0.7985  4.4307 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 2.793    1.671   
##  Residual             2.364    1.537   
## Number of obs: 9900, groups:  id, 28
## 
## Fixed effects:
##                 Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)    6.340e-01  1.568e+00 2.501e+01   0.404   0.6894    
## ImageSetmeth   2.169e+00  3.785e-02 9.870e+03  57.310   <2e-16 ***
## ImageSetopioid 2.165e+00  3.785e-02 9.870e+03  57.198   <2e-16 ***
## Age            1.560e-02  3.968e-02 2.500e+01   0.393   0.6975    
## SexMale        2.218e+00  6.390e-01 2.500e+01   3.471   0.0019 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ImgStm ImgStp Age   
## ImageSetmth -0.012                     
## ImageSetopd -0.012  0.500              
## Age         -0.951  0.000  0.000       
## SexMale     -0.219  0.000  0.000 -0.014
## Type III Analysis of Variance Table with Satterthwaite's method
##           Sum Sq Mean Sq NumDF DenDF   F value    Pr(>F)    
## ImageSet 10330.7  5165.3     2  9870 2185.3590 < 2.2e-16 ***
## Age          0.4     0.4     1    25    0.1546  0.697523    
## Sex         28.5    28.5     1    25   12.0468  0.001899 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = Value ~ ImageSet + Age + Sex + (1 | id), data = dset[dset$RatingType == 
##     RatingType, ])
## 
## Linear Hypotheses:
##                        Estimate Std. Error z value Pr(>|z|)    
## meth - control == 0    2.169091   0.037848  57.310   <1e-06 ***
## opioid - control == 0  2.164848   0.037848  57.198   <1e-06 ***
## opioid - meth == 0    -0.004242   0.037848  -0.112    0.993    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
summarize_one_lme_allimages(subject_data, 'valence')
## [1] "LME OF ALL IMAGES FOR: valence"

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ ImageSet + Age + Sex + (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 34641.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2479 -0.4175  0.0989  0.4949  4.1855 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 1.038    1.019   
##  Residual             1.908    1.381   
## Number of obs: 9897, groups:  id, 28
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       3.12136    0.95728   25.02159   3.261   0.0032 ** 
## ImageSetmeth     -0.44533    0.03402 9867.00911 -13.092   <2e-16 ***
## ImageSetopioid   -0.44667    0.03401 9867.00900 -13.134   <2e-16 ***
## Age               0.03840    0.02423   25.00064   1.585   0.1256    
## SexMale           0.53025    0.39011   25.00638   1.359   0.1862    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ImgStm ImgStp Age   
## ImageSetmth -0.018                     
## ImageSetopd -0.018  0.500              
## Age         -0.951  0.000  0.000       
## SexMale     -0.219  0.000  0.000 -0.014
## Type III Analysis of Variance Table with Satterthwaite's method
##          Sum Sq Mean Sq NumDF DenDF  F value Pr(>F)    
## ImageSet 437.55 218.773     2  9867 114.6380 <2e-16 ***
## Age        4.79   4.794     1    25   2.5120 0.1256    
## Sex        3.53   3.526     1    25   1.8475 0.1862    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = Value ~ ImageSet + Age + Sex + (1 | id), data = dset[dset$RatingType == 
##     RatingType, ])
## 
## Linear Hypotheses:
##                        Estimate Std. Error z value Pr(>|z|)    
## meth - control == 0   -0.445328   0.034017 -13.092   <1e-07 ***
## opioid - control == 0 -0.446667   0.034009 -13.134   <1e-07 ***
## opioid - meth == 0    -0.001339   0.034017  -0.039    0.999    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
summarize_one_lme_allimages(subject_data, 'valence_fromneutral')
## [1] "LME OF ALL IMAGES FOR: valence_fromneutral"

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ ImageSet + Age + Sex + (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 30388.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9879 -0.7091 -0.2381  0.4710  4.0593 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.7221   0.8498  
##  Residual             1.2413   1.1141  
## Number of obs: 9897, groups:  id, 28
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)      -0.39794    0.79833   25.01834  -0.498    0.623    
## ImageSetmeth      0.83280    0.02743 9867.00669  30.357   <2e-16 ***
## ImageSetopioid    0.78788    0.02743 9867.00658  28.726   <2e-16 ***
## Age               0.02146    0.02020   24.99875   1.062    0.298    
## SexMale           0.20780    0.32534   25.00413   0.639    0.529    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ImgStm ImgStp Age   
## ImageSetmth -0.017                     
## ImageSetopd -0.017  0.500              
## Age         -0.951  0.000  0.000       
## SexMale     -0.219  0.000  0.000 -0.014
## Type III Analysis of Variance Table with Satterthwaite's method
##           Sum Sq Mean Sq NumDF DenDF  F value Pr(>F)    
## ImageSet 1447.71  723.86     2  9867 583.1594 <2e-16 ***
## Age         1.40    1.40     1    25   1.1277 0.2984    
## Sex         0.51    0.51     1    25   0.4080 0.5288    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = Value ~ ImageSet + Age + Sex + (1 | id), data = dset[dset$RatingType == 
##     RatingType, ])
## 
## Linear Hypotheses:
##                       Estimate Std. Error z value Pr(>|z|)    
## meth - control == 0    0.83280    0.02743  30.357   <1e-04 ***
## opioid - control == 0  0.78788    0.02743  28.726   <1e-04 ***
## opioid - meth == 0    -0.04493    0.02743  -1.638     0.23    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

new LMEs for 7/2019

drug images: craving|methToOpioid ~ cat drug + time * session

Need category to have the same values for opioid/meth here–e.g. instrument_hand instead of opioid_instrument_hand

#meth_instrument opioid_injection_instrument
#meth_hand opioid_hand
#meth_injection_hand opioid_injection_hand
#meth_instrument_hand opioid_instrument_hand
#meth_face_activities opioid_face_activities
#meth opioid

subject_data$UnifiedCategory <- NA
subject_data$UnifiedCategory[subject_data$Category %in% c('meth', 'opioid')] <- 'drug'
subject_data$UnifiedCategory[subject_data$Category %in% c('meth_instrument', 'opioid_injection_instrument')] <- 'instrument'
subject_data$UnifiedCategory[subject_data$Category %in% c('meth_hand', 'opioid_hand')] <- 'drug_hand'
subject_data$UnifiedCategory[subject_data$Category %in% c('meth_injection_hand', 'opioid_injection_hand')] <- 'injection_hand'
subject_data$UnifiedCategory[subject_data$Category %in% c('meth_instrument_hand', 'opioid_instrument_hand')] <- 'instrument_hand'
subject_data$UnifiedCategory[subject_data$Category %in% c('meth_face_activities', 'opioid_face_activities')] <- 'face_activities'


subject_data$MethToOpioid <- NA
subject_data$MethToOpioid[subject_data$RelatedCategory == 'Meth'] <- -1
subject_data$MethToOpioid[subject_data$RelatedCategory == 'Neither'] <- 0
subject_data$MethToOpioid[subject_data$RelatedCategory == 'Opioid'] <- 1
subject_data$MethToOpioid[subject_data$RelatedCategory == 'Both'] <- 0


summarize_one_lme_drugcats <- function(dset, RatingType){
  this_lme <- lmer(Value ~  UnifiedCategory * ImageSet + Time * Visit + (1 | id), 
                data = dset[dset$RatingType == RatingType,])
  print(paste("LME OF Drug Images FOR:", RatingType))
  print(plot_model(this_lme, title = paste(RatingType), show.values = TRUE, type = 'std'))
  print(summary(this_lme))
  print(anova(this_lme))
}
summarize_one_lme_drugcats(subject_data[subject_data$ImageSet %in% c('meth', 'opioid'),], 'craving')
## [1] "LME OF Drug Images FOR: craving"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: Value ~ UnifiedCategory * ImageSet + Time * Visit + (1 | id)
##    Data: dset[dset$RatingType == RatingType, ]
## 
## REML criterion at convergence: 53644.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.4418 -0.1264  0.1650  0.4942  3.4449 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 132.5    11.51   
##  Residual             194.5    13.95   
## Number of obs: 6600, groups:  id, 28
## 
## Fixed effects:
##                                                 Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                                    9.420e+01  2.299e+00  3.329e+01  40.973  < 2e-16 ***
## UnifiedCategorydrug_hand                      -1.988e+00  9.711e-01  6.558e+03  -2.048 0.040634 *  
## UnifiedCategoryface_activities                -1.302e+00  8.412e-01  6.558e+03  -1.548 0.121746    
## UnifiedCategoryinjection_hand                 -6.386e-01  9.714e-01  6.558e+03  -0.657 0.510990    
## UnifiedCategoryinstrument                     -2.950e+00  7.422e-01  6.558e+03  -3.975 7.13e-05 ***
## UnifiedCategoryinstrument_hand                -2.969e+00  8.415e-01  6.558e+03  -3.528 0.000422 ***
## ImageSetopioid                                -3.334e+00  8.411e-01  6.558e+03  -3.964 7.45e-05 ***
## Time                                          -2.831e-02  4.519e-03  6.558e+03  -6.264 4.00e-10 ***
## VisitV2                                       -1.484e+00  6.910e-01  6.559e+03  -2.148 0.031770 *  
## UnifiedCategorydrug_hand:ImageSetopioid        1.079e+00  1.373e+00  6.558e+03   0.786 0.431959    
## UnifiedCategoryface_activities:ImageSetopioid  4.394e+00  1.189e+00  6.558e+03   3.694 0.000222 ***
## UnifiedCategoryinjection_hand:ImageSetopioid   4.739e+00  1.374e+00  6.558e+03   3.450 0.000564 ***
## UnifiedCategoryinstrument:ImageSetopioid       4.280e+00  1.049e+00  6.558e+03   4.081 4.55e-05 ***
## UnifiedCategoryinstrument_hand:ImageSetopioid  6.292e+00  1.190e+00  6.558e+03   5.289 1.27e-07 ***
## Time:VisitV2                                  -2.774e-02  6.474e-03  6.558e+03  -4.284 1.86e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(summary(this_lme), correlation=TRUE)  or
##     vcov(summary(this_lme))        if you need it

## Type III Analysis of Variance Table with Satterthwaite's method
##                          Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## UnifiedCategory            5567    1113     5 6558.0   5.7243 2.821e-05 ***
## ImageSet                     24      24     1 6558.0   0.1245   0.72422    
## Time                      32938   32938     1 6558.2 169.3537 < 2.2e-16 ***
## Visit                       897     897     1 6559.0   4.6128   0.03177 *  
## UnifiedCategory:ImageSet   7214    1443     5 6558.0   7.4185 5.999e-07 ***
## Time:Visit                 3569    3569     1 6558.1  18.3516 1.863e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dset <- subject_data[subject_data$ImageSet %in% c('meth', 'opioid'),]
this_lme <- lmer(MethToOpioid ~  UnifiedCategory * ImageSet + Time * Visit + (1 | id), 
              data = dset)
print(paste("LME OF Drug Images FOR:", 'MethToOpioid'))
## [1] "LME OF Drug Images FOR: MethToOpioid"
print(plot_model(this_lme, title = paste('MethToOpioid'), show.values = TRUE, type = 'std'))

print(summary(this_lme))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MethToOpioid ~ UnifiedCategory * ImageSet + Time * Visit + (1 |      id)
##    Data: dset
## 
## REML criterion at convergence: 11106.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3056 -0.6445 -0.2287  0.8722  3.7004 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.02873  0.1695  
##  Residual             0.30705  0.5541  
## Number of obs: 6597, groups:  id, 28
## 
## Fixed effects:
##                                                 Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                                   -1.033e+00  4.361e-02  8.478e+01 -23.693  < 2e-16 ***
## UnifiedCategorydrug_hand                       3.599e-01  3.862e-02  6.555e+03   9.319  < 2e-16 ***
## UnifiedCategoryface_activities                 3.032e-01  3.342e-02  6.555e+03   9.072  < 2e-16 ***
## UnifiedCategoryinjection_hand                  8.629e-01  3.864e-02  6.555e+03  22.336  < 2e-16 ***
## UnifiedCategoryinstrument                      3.153e-01  2.949e-02  6.555e+03  10.690  < 2e-16 ***
## UnifiedCategoryinstrument_hand                 1.660e-01  3.343e-02  6.555e+03   4.965 7.05e-07 ***
## ImageSetopioid                                 1.585e+00  3.342e-02  6.555e+03  47.415  < 2e-16 ***
## Time                                           4.816e-04  1.796e-04  6.556e+03   2.682  0.00733 ** 
## VisitV2                                        1.357e-01  2.745e-02  6.561e+03   4.943 7.89e-07 ***
## UnifiedCategorydrug_hand:ImageSetopioid       -5.983e-01  5.459e-02  6.555e+03 -10.959  < 2e-16 ***
## UnifiedCategoryface_activities:ImageSetopioid -7.466e-01  4.726e-02  6.555e+03 -15.796  < 2e-16 ***
## UnifiedCategoryinjection_hand:ImageSetopioid  -1.336e+00  5.460e-02  6.555e+03 -24.466  < 2e-16 ***
## UnifiedCategoryinstrument:ImageSetopioid      -6.712e-01  4.168e-02  6.555e+03 -16.101  < 2e-16 ***
## UnifiedCategoryinstrument_hand:ImageSetopioid -5.402e-01  4.726e-02  6.555e+03 -11.430  < 2e-16 ***
## Time:VisitV2                                  -3.814e-04  2.573e-04  6.556e+03  -1.483  0.13821    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(summary(this_lme), correlation=TRUE)  or
##     vcov(summary(this_lme))        if you need it
print(anova(this_lme))
## Type III Analysis of Variance Table with Satterthwaite's method
##                           Sum Sq Mean Sq NumDF  DenDF   F value    Pr(>F)    
## UnifiedCategory            45.01    9.00     5 6555.1   29.3151 < 2.2e-16 ***
## ImageSet                 1257.39 1257.39     1 6555.1 4095.0401 < 2.2e-16 ***
## Time                        1.57    1.57     1 6556.1    5.1050   0.02389 *  
## Visit                       7.50    7.50     1 6561.4   24.4316 7.891e-07 ***
## UnifiedCategory:ImageSet  197.43   39.49     5 6555.1  128.5940 < 2.2e-16 ***
## Time:Visit                  0.67    0.67     1 6556.0    2.1983   0.13821    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

RT to craving ratings

summarize_one_lme_time <- function(dset, formula_to_use){
  this_lme <- lmer(formula_to_use, 
                data = dset)
  print(paste("LME OF ALL IMAGES FOR: Time On Page"))
  print(plot_model(this_lme, title = paste('Time On Page for Craving Ratings'), show.values = TRUE, type = 'std', order.terms = c(3,4,1,2)))
  print(summary(this_lme))
  print(anova(this_lme))
  
  this_lme.ht <- glht(this_lme, linfct = mcp(ImageSet = "Tukey"))
  summary(this_lme.ht, test=adjusted("none"))
  summary(this_lme.ht)
}

craving_data <- subject_data[subject_data$RatingType == 'craving',]

data_touse <- craving_data[craving_data$millisecondsOnPage < 20000,]
print(paste0("Fraction Excluded for >20s on page: ", 1 - (nrow(data_touse) / nrow(craving_data))))
## [1] "Fraction Excluded for >20s on page: 0.0060606060606061"
summarize_one_lme_time(data_touse, as.formula('millisecondsOnPage / 1000 ~  ImageSet + Age + Sex +(1 | id)'))
## [1] "LME OF ALL IMAGES FOR: Time On Page"

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_to_use
##    Data: dset
## 
## REML criterion at convergence: 39808.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1038 -0.5020 -0.2116  0.1848  9.5495 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.7481   0.865   
##  Residual             3.3002   1.817   
## Number of obs: 9840, groups:  id, 28
## 
## Fixed effects:
##                 Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)    1.462e+00  8.159e-01 2.504e+01   1.791 0.085348 .  
## ImageSetmeth   1.745e-01  4.485e-02 9.810e+03   3.890 0.000101 ***
## ImageSetopioid 3.084e-01  4.486e-02 9.810e+03   6.874 6.62e-12 ***
## Age            3.856e-02  2.064e-02 2.499e+01   1.868 0.073535 .  
## SexMale        3.581e-01  3.324e-01 2.501e+01   1.077 0.291712    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ImgStm ImgStp Age   
## ImageSetmth -0.027                     
## ImageSetopd -0.027  0.500              
## Age         -0.951  0.000  0.000       
## SexMale     -0.219  0.000  0.000 -0.014
## Type III Analysis of Variance Table with Satterthwaite's method
##           Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## ImageSet 156.853  78.426     2  9810 23.7642 5.061e-11 ***
## Age       11.515  11.515     1    25  3.4893   0.07353 .  
## Sex        3.829   3.829     1    25  1.1602   0.29171    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = formula_to_use, data = dset)
## 
## Linear Hypotheses:
##                       Estimate Std. Error z value Pr(>|z|)    
## meth - control == 0    0.17449    0.04485   3.890  < 0.001 ***
## opioid - control == 0  0.30836    0.04486   6.874  < 0.001 ***
## opioid - meth == 0     0.13387    0.04487   2.984  0.00803 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
summarize_one_lme_time(data_touse, as.formula('millisecondsOnPage / 1000 ~  ImageSet + Age + Sex + Time * Visit + (1 | id)'))
## [1] "LME OF ALL IMAGES FOR: Time On Page"
## Number of values in `order.terms` does not match number of terms. Terms are not sorted.

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_to_use
##    Data: dset
## 
## REML criterion at convergence: 38644
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2602 -0.5112 -0.1761  0.2246 10.2432 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.7705   0.8778  
##  Residual             2.9220   1.7094  
## Number of obs: 9840, groups:  id, 28
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     2.825e+00  8.284e-01  2.517e+01   3.410   0.0022 ** 
## ImageSetmeth    2.529e-01  4.231e-02  9.807e+03   5.979 2.33e-09 ***
## ImageSetopioid  3.763e-01  4.229e-02  9.807e+03   8.897  < 2e-16 ***
## Age             3.909e-02  2.093e-02  2.498e+01   1.868   0.0736 .  
## SexMale         3.364e-01  3.371e-01  2.499e+01   0.998   0.3279    
## Time           -1.180e-02  4.664e-04  9.807e+03 -25.305  < 2e-16 ***
## VisitV2        -1.240e+00  6.900e-02  9.809e+03 -17.967  < 2e-16 ***
## Time:VisitV2    5.510e-03  6.642e-04  9.807e+03   8.296  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ImgStm ImgStp Age    SexMal Time   VistV2
## ImageSetmth -0.023                                          
## ImageSetopd -0.024  0.502                                   
## Age         -0.950  0.000  0.000                            
## SexMale     -0.219  0.000  0.000 -0.014                     
## Time        -0.049 -0.046 -0.036  0.000  0.000              
## VisitV2     -0.041  0.005  0.010 -0.001  0.001  0.605       
## Time:VistV2  0.036 -0.005 -0.012  0.000  0.000 -0.700 -0.864
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)    
## ImageSet    240.22  120.11     2 9807.0  41.105 < 2e-16 ***
## Age          10.19   10.19     1   25.0   3.488 0.07359 .  
## Sex           2.91    2.91     1   25.0   0.996 0.32785    
## Time       2155.99 2155.99     1 9807.0 737.839 < 2e-16 ***
## Visit       943.31  943.31     1 9808.8 322.827 < 2e-16 ***
## Time:Visit  201.13  201.13     1 9807.0  68.831 < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = formula_to_use, data = dset)
## 
## Linear Hypotheses:
##                       Estimate Std. Error z value Pr(>|z|)    
## meth - control == 0    0.25294    0.04231   5.979  < 0.001 ***
## opioid - control == 0  0.37628    0.04229   8.897  < 0.001 ***
## opioid - meth == 0     0.12333    0.04222   2.921  0.00976 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Total time to complete task

total_values <- aggregate(subject_data$millisecondsOnPage, 
                            by = list(id = subject_data$id, Visit = subject_data$Visit), 
                            FUN = sum, na.rm = TRUE)

names(total_values) <- c('id', 'Visit', 'TaskCompletionTime')
#convert to seconds
total_values$TaskCompletionTime <- total_values$TaskCompletionTime / 1000

print("Total task completion time for each visit:")
## [1] "Total task completion time for each visit:"
CreateTableOne(vars = 'TaskCompletionTime', data = total_values, strata = 'Visit')
##                                 Stratified by Visit
##                                  V1                V2                p      test
##   n                                   28                27                      
##   TaskCompletionTime (mean (sd)) 4793.03 (1653.44) 3583.79 (1081.61)  0.002

Time On Image as a function of craving rating, just among drug images

image_total_values <- aggregate(subject_data$millisecondsOnPage, 
                            by = list(id = subject_data$id, ImageSetFile = subject_data$ImageSetFile), 
                            FUN = sum, na.rm = TRUE)

names(image_total_values) <- c('id', 'ImageSetFile', 'ImageCompletionTime')
#convert to seconds
image_total_values$ImageCompletionTime <- image_total_values$ImageCompletionTime / 1000

subject_data <- merge(subject_data, image_total_values)


summarize_one_lme_imagetime <- function(dset, formula_to_use){
  this_lme <- lmer(formula_to_use, 
                data = dset)
  print(paste("LME OF DRUG IMAGES FOR: Time On Image"))
  print(plot_model(this_lme, title = paste('Time On Image for Craving Ratings'), show.values = TRUE, type = 'std', order.terms = c(3,4,2,1)))
  print(summary(this_lme))
  print(anova(this_lme))
  
  #this_lme.ht <- glht(this_lme, linfct = mcp(ImageSet = "Tukey"))
  #summary(this_lme.ht, test=adjusted("none"))
  #summary(this_lme.ht)
}

craving_data <- subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet %in% c('meth', 'opioid'),]
names(craving_data)[names(craving_data) == 'Value'] <- 'craving'
data_touse <- craving_data[craving_data$ImageCompletionTime < 100,]

print(paste0("Fraction Excluded for >100s on image: ", 1 - (nrow(data_touse) / nrow(craving_data))))
## [1] "Fraction Excluded for >100s on image: 0.0106060606060606"
summarize_one_lme_imagetime(data_touse, as.formula('ImageCompletionTime ~  craving + ImageSet + Time * Visit + (1 | id)'))
## [1] "LME OF DRUG IMAGES FOR: Time On Image"
## Number of values in `order.terms` does not match number of terms. Terms are not sorted.

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_to_use
##    Data: dset
## 
## REML criterion at convergence: 48833.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2938 -0.5597 -0.1788  0.2675  8.5762 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept)  36.87    6.072  
##  Residual             101.32   10.066  
## Number of obs: 6530, groups:  id, 28
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     2.965e+01  1.459e+00  6.835e+01  20.315  < 2e-16 ***
## craving         3.829e-02  8.876e-03  6.506e+03   4.314 1.63e-05 ***
## ImageSetopioid  1.134e+00  2.492e-01  6.497e+03   4.552 5.40e-06 ***
## Time           -9.147e-02  3.308e-03  6.498e+03 -27.651  < 2e-16 ***
## VisitV2        -7.773e+00  5.042e-01  6.499e+03 -15.415  < 2e-16 ***
## Time:VisitV2    3.021e-02  4.715e-03  6.498e+03   6.408 1.58e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cravng ImgStp Time   VistV2
## craving     -0.562                            
## ImageSetopd -0.082 -0.010                     
## Time        -0.252  0.077  0.010              
## VisitV2     -0.184  0.026  0.007  0.607       
## Time:VistV2  0.117  0.053 -0.009 -0.693 -0.864
## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## craving      1885    1885     1 6506.3   18.608 1.630e-05 ***
## ImageSet     2100    2100     1 6497.0   20.725 5.400e-06 ***
## Time       103901  103901     1 6499.5 1025.490 < 2.2e-16 ***
## Visit       24075   24075     1 6499.0  237.618 < 2.2e-16 ***
## Time:Visit   4160    4160     1 6497.6   41.061 1.579e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_imagetime(data_touse, as.formula('ImageCompletionTime ~  craving + Time * Visit + (1 | id)'))
## [1] "LME OF DRUG IMAGES FOR: Time On Image"

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_to_use
##    Data: dset
## 
## REML criterion at convergence: 48853.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3473 -0.5595 -0.1793  0.2651  8.5086 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept)  36.86    6.071  
##  Residual             101.63   10.081  
## Number of obs: 6530, groups:  id, 28
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   3.019e+01  1.455e+00  6.761e+01  20.746  < 2e-16 ***
## craving       3.870e-02  8.889e-03  6.507e+03   4.353 1.36e-05 ***
## Time         -9.162e-02  3.313e-03  6.499e+03 -27.656  < 2e-16 ***
## VisitV2      -7.789e+00  5.050e-01  6.500e+03 -15.423  < 2e-16 ***
## Time:VisitV2  3.041e-02  4.722e-03  6.499e+03   6.439 1.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cravng Time   VistV2
## craving     -0.566                     
## Time        -0.253  0.077              
## VisitV2     -0.184  0.026  0.607       
## Time:VistV2  0.117  0.053 -0.693 -0.864
## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## craving      1926    1926     1 6507.0   18.951 1.362e-05 ***
## Time       104049  104049     1 6500.5 1023.842 < 2.2e-16 ***
## Visit       24173   24173     1 6500.1  237.863 < 2.2e-16 ***
## Time:Visit   4213    4213     1 6498.6   41.459 1.290e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summarize_one_lme_imagetime(data_touse, as.formula('ImageCompletionTime ~  craving + (1 | id)'))
## [1] "LME OF DRUG IMAGES FOR: Time On Image"
## Number of values in `order.terms` does not match number of terms. Terms are not sorted.

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: formula_to_use
##    Data: dset
## 
## REML criterion at convergence: 50129.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0339 -0.5448 -0.2364  0.2137  7.6031 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept)  35.35    5.946  
##  Residual             124.01   11.136  
## Number of obs: 6530, groups:  id, 28
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) 1.353e+01  1.411e+00 6.443e+01   9.589 5.03e-14 ***
## craving     1.043e-01  9.574e-03 6.489e+03  10.897  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## craving -0.597
## Type III Analysis of Variance Table with Satterthwaite's method
##         Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## craving  14724   14724     1 6489.4  118.73 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

RT LMEs with clinical variables

this_lme <- lmer(millisecondsOnPage / 1000 ~ Age + Sex + MethOnsetAge + MethMonthsUsed + MethDollars + MethDaysLastMonth + MethDaysAbstinant + Time * Visit +(1 | id), 
              data = subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet == 'meth',])
## Warning: Some predictor variables are on very different scales: consider rescaling

## Warning: Some predictor variables are on very different scales: consider rescaling
print(plot_model(this_lme, title = 'meth RT to Craving Ratings', show.values = TRUE, type = 'std', order.terms = c(1,2,8,9,10,3,4,6,5,7)))

print(summary(this_lme))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: millisecondsOnPage/1000 ~ Age + Sex + MethOnsetAge + MethMonthsUsed +  
##     MethDollars + MethDaysLastMonth + MethDaysAbstinant + Time *      Visit + (1 | id)
##    Data: subject_data[subject_data$RatingType == "craving" & subject_data$ImageSet ==  
##     "meth", ]
## 
## REML criterion at convergence: 17740.6
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -1.012 -0.293 -0.116  0.071 32.336 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept)  0.7159  0.8461  
##  Residual             14.9185  3.8625  
## Number of obs: 3180, groups:  id, 27
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        4.341e+00  1.529e+00  1.980e+01   2.840   0.0102 *  
## Age                5.741e-02  2.813e-02  1.918e+01   2.040   0.0553 .  
## SexMale           -3.078e-01  4.682e-01  1.916e+01  -0.657   0.5188    
## MethOnsetAge       2.185e-02  4.571e-02  1.922e+01   0.478   0.6381    
## MethMonthsUsed    -6.900e-04  2.857e-03  1.905e+01  -0.241   0.8118    
## MethDollars        2.532e-06  9.471e-05  1.946e+01   0.027   0.9789    
## MethDaysLastMonth -3.776e-02  1.845e-02  1.925e+01  -2.047   0.0546 .  
## MethDaysAbstinant -1.063e-03  6.759e-04  1.908e+01  -1.573   0.1323    
## Time              -1.366e-02  1.804e-03  3.153e+03  -7.569 4.89e-14 ***
## VisitV2           -1.305e+00  2.778e-01  3.158e+03  -4.697 2.76e-06 ***
## Time:VisitV2       4.547e-03  2.591e-03  3.153e+03   1.755   0.0793 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Age    SexMal MthOnA MthMnU MthDll MthDLM MthDyA Time   VistV2
## Age         -0.526                                                               
## SexMale     -0.053 -0.007                                                        
## MethOnsetAg -0.570 -0.147 -0.390                                                 
## MthMnthsUsd -0.216 -0.386 -0.325  0.573                                          
## MethDollars -0.341  0.121 -0.017  0.174 -0.058                                   
## MthDysLstMn -0.145 -0.115  0.350 -0.023 -0.069 -0.438                            
## MthDysAbstn -0.504  0.076  0.329  0.055  0.098  0.390  0.053                     
## Time        -0.107 -0.001 -0.007  0.006  0.007 -0.006 -0.003 -0.013              
## VisitV2     -0.075 -0.012  0.006 -0.007  0.010 -0.010  0.003 -0.010  0.602       
## Time:VistV2  0.071  0.006  0.000 -0.001 -0.006 -0.001  0.005  0.005 -0.695 -0.868
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
print(anova(this_lme))
## Type III Analysis of Variance Table with Satterthwaite's method
##                    Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## Age                 62.12   62.12     1   19.18  4.1636   0.05531 .  
## Sex                  6.45    6.45     1   19.16  0.4321   0.51880    
## MethOnsetAge         3.41    3.41     1   19.22  0.2284   0.63808    
## MethMonthsUsed       0.87    0.87     1   19.05  0.0583   0.81175    
## MethDollars          0.01    0.01     1   19.46  0.0007   0.97894    
## MethDaysLastMonth   62.50   62.50     1   19.25  4.1893   0.05458 .  
## MethDaysAbstinant   36.89   36.89     1   19.08  2.4730   0.13225    
## Time              1148.29 1148.29     1 3154.39 76.9707 < 2.2e-16 ***
## Visit              329.06  329.06     1 3158.00 22.0574 2.759e-06 ***
## Time:Visit          45.96   45.96     1 3152.70  3.0805   0.07934 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
this_lme <- lmer(millisecondsOnPage / 1000 ~ Age + Sex + OpioidOrHeroinOnsetAge + OpioidOrHeroinMonthsUsed + OpioidOrHeroinDollars + 
                   OpioidOrHeroinDaysLastMonth + OpioidOrHeroinDaysAbstinant + Time * Visit + (1 | id), 
              data = subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet == 'opioid',])
## Warning: Some predictor variables are on very different scales: consider rescaling

## Warning: Some predictor variables are on very different scales: consider rescaling
print(plot_model(this_lme, title = 'opioid RT to Craving Ratings', show.values = TRUE, type = 'std', order.terms = c(1,2,8, 9, 10, 3,4,6,5,7)))

print(summary(this_lme))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: millisecondsOnPage/1000 ~ Age + Sex + OpioidOrHeroinOnsetAge +  
##     OpioidOrHeroinMonthsUsed + OpioidOrHeroinDollars + OpioidOrHeroinDaysLastMonth +  
##     OpioidOrHeroinDaysAbstinant + Time * Visit + (1 | id)
##    Data: subject_data[subject_data$RatingType == "craving" & subject_data$ImageSet ==  
##     "opioid", ]
## 
## REML criterion at convergence: 17831.3
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -1.026 -0.323 -0.127  0.084 36.956 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept)  1.056   1.028   
##  Residual             15.314   3.913   
## Number of obs: 3180, groups:  id, 27
## 
## Fixed effects:
##                               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                  4.250e+00  1.500e+00  1.939e+01   2.833   0.0105 *  
## Age                          4.577e-02  4.288e-02  1.891e+01   1.067   0.2993    
## SexMale                      2.286e-01  4.654e-01  1.912e+01   0.491   0.6289    
## OpioidOrHeroinOnsetAge       1.486e-02  4.888e-02  1.920e+01   0.304   0.7644    
## OpioidOrHeroinMonthsUsed    -6.160e-04  3.490e-03  1.903e+01  -0.176   0.8618    
## OpioidOrHeroinDollars       -1.590e-04  1.643e-04  1.918e+01  -0.968   0.3451    
## OpioidOrHeroinDaysLastMonth -1.926e-02  2.604e-02  1.887e+01  -0.740   0.4686    
## OpioidOrHeroinDaysAbstinant -1.380e-04  1.879e-04  1.894e+01  -0.735   0.4715    
## Time                        -1.509e-02  1.833e-03  3.153e+03  -8.235 2.61e-16 ***
## VisitV2                     -1.478e+00  2.772e-01  3.157e+03  -5.332 1.04e-07 ***
## Time:VisitV2                 5.449e-03  2.611e-03  3.152e+03   2.087   0.0370 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Age    SexMal OpOHOA OpOHMU OpdOHD OOHDLM OpOHDA Time   VistV2
## Age         -0.750                                                               
## SexMale     -0.065  0.021                                                        
## OpdOrHrnOnA  0.083 -0.616 -0.149                                                 
## OpdOrHrnMnU  0.125 -0.397 -0.336  0.409                                          
## OpdOrHrnDll  0.058  0.119  0.215 -0.421 -0.250                                   
## OpdOrHrnDLM -0.629  0.486  0.056 -0.261 -0.392 -0.157                            
## OpdOrHrnDyA -0.177 -0.045 -0.160 -0.003  0.227  0.016  0.236                     
## Time        -0.104 -0.002  0.003 -0.006 -0.003  0.004 -0.005  0.001              
## VisitV2     -0.080 -0.004  0.006 -0.007 -0.008  0.010 -0.005 -0.003  0.603       
## Time:VistV2  0.069  0.008  0.002 -0.002  0.002 -0.002  0.005 -0.002 -0.702 -0.863
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
print(anova(this_lme))
## Type III Analysis of Variance Table with Satterthwaite's method
##                              Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## Age                           17.44   17.44     1   18.91  1.1390   0.29930    
## Sex                            3.70    3.70     1   19.12  0.2413   0.62886    
## OpioidOrHeroinOnsetAge         1.41    1.41     1   19.20  0.0924   0.76444    
## OpioidOrHeroinMonthsUsed       0.48    0.48     1   19.03  0.0311   0.86177    
## OpioidOrHeroinDollars         14.35   14.35     1   19.18  0.9371   0.34509    
## OpioidOrHeroinDaysLastMonth    8.38    8.38     1   18.87  0.5471   0.46861    
## OpioidOrHeroinDaysAbstinant    8.27    8.27     1   18.94  0.5399   0.47146    
## Time                        1374.11 1374.11     1 3152.11 89.7303 < 2.2e-16 ***
## Visit                        435.33  435.33     1 3156.68 28.4271 1.041e-07 ***
## Time:Visit                    66.69   66.69     1 3152.34  4.3547   0.03699 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
this_lme <- lmer( ImageCompletionTime ~ Age + Sex + MethOnsetAge + MethMonthsUsed + MethDollars + MethDaysLastMonth + MethDaysAbstinant + Time * Visit +(1 | id), 
              data = subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet == 'meth',])
## Warning: Some predictor variables are on very different scales: consider rescaling

## Warning: Some predictor variables are on very different scales: consider rescaling
print(plot_model(this_lme, title = 'meth Image Completion Time', show.values = TRUE, type = 'std', order.terms = c(1,2,8,9,10,3,4,6,5,7)))

print(summary(this_lme))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: ImageCompletionTime ~ Age + Sex + MethOnsetAge + MethMonthsUsed +  
##     MethDollars + MethDaysLastMonth + MethDaysAbstinant + Time *      Visit + (1 | id)
##    Data: subject_data[subject_data$RatingType == "craving" & subject_data$ImageSet ==  
##     "meth", ]
## 
## REML criterion at convergence: 30424.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.8646 -0.2624 -0.1015  0.0871 31.3912 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept)  48.59    6.971  
##  Residual             815.56   28.558  
## Number of obs: 3180, groups:  id, 27
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        2.277e+01  1.239e+01  1.957e+01   1.838   0.0813 .  
## Age                3.293e-01  2.283e-01  1.906e+01   1.442   0.1655    
## SexMale           -2.542e+00  3.800e+00  1.904e+01  -0.669   0.5116    
## MethOnsetAge       3.369e-01  3.709e-01  1.909e+01   0.908   0.3751    
## MethMonthsUsed     7.553e-03  2.320e-02  1.895e+01   0.326   0.7483    
## MethDollars       -8.147e-05  7.681e-04  1.930e+01  -0.106   0.9166    
## MethDaysLastMonth -3.635e-02  1.497e-01  1.911e+01  -0.243   0.8108    
## MethDaysAbstinant -3.952e-03  5.487e-03  1.897e+01  -0.720   0.4801    
## Time              -1.059e-01  1.334e-02  3.153e+03  -7.943 2.73e-15 ***
## VisitV2           -1.116e+01  2.054e+00  3.157e+03  -5.435 5.91e-08 ***
## Time:VisitV2       3.938e-02  1.916e-02  3.152e+03   2.056   0.0399 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Age    SexMal MthOnA MthMnU MthDll MthDLM MthDyA Time   VistV2
## Age         -0.526                                                               
## SexMale     -0.054 -0.006                                                        
## MethOnsetAg -0.570 -0.148 -0.389                                                 
## MthMnthsUsd -0.217 -0.386 -0.326  0.575                                          
## MethDollars -0.340  0.120 -0.015  0.172 -0.057                                   
## MthDysLstMn -0.147 -0.114  0.349 -0.022 -0.070 -0.438                            
## MthDysAbstn -0.505  0.075  0.330  0.054  0.099  0.389  0.054                     
## Time        -0.098 -0.001 -0.006  0.005  0.006 -0.006 -0.003 -0.012              
## VisitV2     -0.068 -0.011  0.006 -0.006  0.009 -0.009  0.003 -0.009  0.602       
## Time:VistV2  0.065  0.005  0.000 -0.001 -0.005 -0.001  0.005  0.005 -0.695 -0.868
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
print(anova(this_lme))
## Type III Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## Age                 1696    1696     1   19.06  2.0798   0.16550    
## Sex                  365     365     1   19.04  0.4474   0.51161    
## MethOnsetAge         673     673     1   19.09  0.8249   0.37508    
## MethMonthsUsed        86      86     1   18.95  0.1060   0.74828    
## MethDollars            9       9     1   19.30  0.0112   0.91663    
## MethDaysLastMonth     48      48     1   19.11  0.0589   0.81076    
## MethDaysAbstinant    423     423     1   18.97  0.5189   0.48010    
## Time               65945   65945     1 3153.62 80.8593 < 2.2e-16 ***
## Visit              24087   24087     1 3156.88 29.5346  5.91e-08 ***
## Time:Visit          3447    3447     1 3152.17  4.2264   0.03988 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
this_lme <- lmer(ImageCompletionTime ~ Age + Sex + OpioidOrHeroinOnsetAge + OpioidOrHeroinMonthsUsed + OpioidOrHeroinDollars + 
                   OpioidOrHeroinDaysLastMonth + OpioidOrHeroinDaysAbstinant + Time * Visit + (1 | id), 
              data = subject_data[subject_data$RatingType == 'craving' & subject_data$ImageSet == 'opioid',])
## Warning: Some predictor variables are on very different scales: consider rescaling

## Warning: Some predictor variables are on very different scales: consider rescaling
print(plot_model(this_lme, title = 'opioid Image Completion Time', show.values = TRUE, type = 'std', order.terms = c(1,2,8, 9, 10, 3,4,6,5,7)))

print(summary(this_lme))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: ImageCompletionTime ~ Age + Sex + OpioidOrHeroinOnsetAge + OpioidOrHeroinMonthsUsed +  
##     OpioidOrHeroinDollars + OpioidOrHeroinDaysLastMonth + OpioidOrHeroinDaysAbstinant +  
##     Time * Visit + (1 | id)
##    Data: subject_data[subject_data$RatingType == "craving" & subject_data$ImageSet ==  
##     "opioid", ]
## 
## REML criterion at convergence: 29457.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2662 -0.3552 -0.1292  0.1051 26.5328 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept)  65.71    8.106  
##  Residual             598.76   24.470  
## Number of obs: 3180, groups:  id, 27
## 
## Fixed effects:
##                               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                  3.465e+01  1.157e+01  1.929e+01   2.996  0.00734 ** 
## Age                          5.175e-02  3.314e-01  1.897e+01   0.156  0.87756    
## SexMale                      9.412e-01  3.592e+00  1.912e+01   0.262  0.79611    
## OpioidOrHeroinOnsetAge       4.239e-01  3.772e-01  1.918e+01   1.124  0.27501    
## OpioidOrHeroinMonthsUsed     9.660e-03  2.695e-02  1.906e+01   0.358  0.72398    
## OpioidOrHeroinDollars       -1.931e-03  1.268e-03  1.916e+01  -1.523  0.14404    
## OpioidOrHeroinDaysLastMonth -6.221e-02  2.012e-01  1.895e+01  -0.309  0.76058    
## OpioidOrHeroinDaysAbstinant -4.398e-04  1.451e-03  1.899e+01  -0.303  0.76519    
## Time                        -1.401e-01  1.146e-02  3.152e+03 -12.226  < 2e-16 ***
## VisitV2                     -1.153e+01  1.734e+00  3.155e+03  -6.648 3.49e-11 ***
## Time:VisitV2                 4.732e-02  1.633e-02  3.152e+03   2.898  0.00378 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Age    SexMal OpOHOA OpOHMU OpdOHD OOHDLM OpOHDA Time   VistV2
## Age         -0.752                                                               
## SexMale     -0.066  0.021                                                        
## OpdOrHrnOnA  0.083 -0.616 -0.146                                                 
## OpdOrHrnMnU  0.126 -0.396 -0.334  0.407                                          
## OpdOrHrnDll  0.058  0.119  0.212 -0.419 -0.248                                   
## OpdOrHrnDLM -0.631  0.486  0.056 -0.262 -0.393 -0.157                            
## OpdOrHrnDyA -0.177 -0.045 -0.159 -0.005  0.227  0.017  0.236                     
## Time        -0.084 -0.002  0.003 -0.005 -0.003  0.003 -0.004  0.001              
## VisitV2     -0.065 -0.003  0.005 -0.006 -0.006  0.008 -0.004 -0.002  0.603       
## Time:VistV2  0.056  0.006  0.002 -0.002  0.001 -0.001  0.004 -0.001 -0.702 -0.863
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
print(anova(this_lme))
## Type III Analysis of Variance Table with Satterthwaite's method
##                             Sum Sq Mean Sq NumDF   DenDF  F value    Pr(>F)    
## Age                             15      15     1   18.97   0.0244  0.877562    
## Sex                             41      41     1   19.12   0.0686  0.796114    
## OpioidOrHeroinOnsetAge         756     756     1   19.18   1.2627  0.275008    
## OpioidOrHeroinMonthsUsed        77      77     1   19.06   0.1285  0.723983    
## OpioidOrHeroinDollars         1389    1389     1   19.16   2.3202  0.144038    
## OpioidOrHeroinDaysLastMonth     57      57     1   18.95   0.0956  0.760579    
## OpioidOrHeroinDaysAbstinant     55      55     1   18.99   0.0918  0.765194    
## Time                        121834  121834     1 3151.43 203.4770 < 2.2e-16 ***
## Visit                        26463   26463     1 3154.74  44.1971 3.485e-11 ***
## Time:Visit                    5029    5029     1 3151.58   8.3997  0.003779 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Raincloud Plots of RTs after excluding some

rc_data1 <- subject_data[subject_data$millisecondsOnPage < 20000 & subject_data$RatingType == 'craving',]
p4 <- ggplot(rc_data1,aes(x=ImageSetCap,y=millisecondsOnPage / 1000, fill = ImageSetCap, colour = ImageSetCap))+
  geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim =  FALSE)+
  geom_point(position = position_jitter(width = .15), size = .25)+
  geom_boxplot(aes(x = as.numeric(ImageSetCap)+0.25, y = millisecondsOnPage / 1000),outlier.shape = NA,
  alpha = 0.3, width = .1, colour = "BLACK") +
  ylab('')+xlab('')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
  colour = FALSE) +
  scale_color_manual(values = cols) +
  scale_fill_manual(values = cols) +
  ggtitle("SecondsOnPage: Craving Ratings") #+ scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9))

rc_data2 <- subject_data[subject_data$ImageCompletionTime < 100,]
p5 <- ggplot(rc_data2,aes(x=ImageSetCap,y=ImageCompletionTime, fill = ImageSetCap, colour = ImageSetCap))+
  geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2, trim =  FALSE)+
  geom_point(position = position_jitter(width = .15), size = .25)+
  geom_boxplot(aes(x = as.numeric(ImageSetCap)+0.25, y = ImageCompletionTime),outlier.shape = NA,
  alpha = 0.3, width = .1, colour = "BLACK") +
  ylab('')+xlab('')+coord_flip()+theme_cowplot()+guides(fill = FALSE,
  colour = FALSE) +
  scale_color_manual(values = cols) +
  scale_fill_manual(values = cols) +
  ggtitle("ImageCompletionTime") #+ scale_y_continuous(breaks = c(1, 3, 5, 7, 9), limits = c(1, 9))

png('RTraincloudsExcluded.png', width = 1400, height = 400)
print(annotate_figure(ggarrange(p4, p5, ncol = 1, nrow = 2, common.legend = TRUE, legend = 'right'), top = text_grob('RTRainclouds', color = 'red'))) 
dev.off()
## png 
##   2

Make random sets of images

#remove group summary lines
selection_data <- merged_summary_values[!is.na(merged_summary_values$File),]

#remove opioid image 118 because 109 and 118 are the same 
selection_data <- selection_data[(selection_data$Group != 'opioid Slide118.jpeg'),]

#remove images with text on them
selection_data <- selection_data[!(selection_data$Group %in%  c('opioid Slide015.jpeg', 'opioid Slide016.jpeg', 'opioid Slide017.jpeg', 
                                                                'opioid Slide018.jpeg', 'opioid Slide019.jpeg', 'opioid Slide020.jpeg',
                                                                'control Slide017.jpeg', 'control Slide018.jpeg', 'control Slide019.jpeg', 'control Slide020.jpeg')),]

#remove light bulb images, since they can be associated with drug use
selection_data <- selection_data[!(selection_data$Group %in%  c('opioid Slide57.jpeg', 'opioid Slide87.jpeg')),]




#line below to select data based on hightst/lowest craving
#sorted_data <- selection_data[order(selection_data$allsubjects_craving_mean),]

#instead, let's select based on hue and see if that helps balance things--get images closest to 0.2
try_one_criteria <- function(control_target_hue, control_target_saturation, control_target_value,
                             meth_target_hue, meth_target_saturation, meth_target_value,
                             opioid_target_hue, opioid_target_saturation, opioid_target_value){
  #sorted_data <- selection_data[sample(nrow(selection_data)),]
  
  selection_data$selection_cost <- NA
  selection_data$selection_cost[selection_data$ImageSet == 'control'] <- 
    (selection_data$hue_mean[selection_data$ImageSet == 'control'] - control_target_hue)^2 + 
    (selection_data$saturation_mean[selection_data$ImageSet == 'control'] - control_target_saturation)^2 +
    (selection_data$value_mean[selection_data$ImageSet == 'control'] - control_target_value)^2
  
  selection_data$selection_cost[selection_data$ImageSet == 'meth'] <- 
    (selection_data$hue_mean[selection_data$ImageSet == 'meth'] - meth_target_hue)^2 + 
    (selection_data$saturation_mean[selection_data$ImageSet == 'meth'] - meth_target_saturation)^2 +
    (selection_data$value_mean[selection_data$ImageSet == 'meth'] - meth_target_value)^2
  
  selection_data$selection_cost[selection_data$ImageSet == 'opioid'] <- 
    (selection_data$hue_mean[selection_data$ImageSet == 'opioid'] - opioid_target_hue)^2 + 
    (selection_data$saturation_mean[selection_data$ImageSet == 'opioid'] - opioid_target_saturation)^2 +
    (selection_data$value_mean[selection_data$ImageSet == 'opioid'] - opioid_target_value)^2
  
  
  (selection_data$hue_mean - 0.3)^2 + (selection_data$saturation_mean - 0.3)^2 + (selection_data$value_mean - 0.3)^2
  sorted_data <- selection_data[order(selection_data$selection_cost),]


  
  #label categories 1:6, to make indexind into number_selected easier
  
  sorted_data$category_number <- NA
  
  sorted_data$category_number[sorted_data$Category == 'control_objects'] <- 1
  sorted_data$category_number[sorted_data$Category == 'control_objects_hand'] <- 2
  sorted_data$category_number[sorted_data$Category == 'control_tool'] <- 3
  sorted_data$category_number[sorted_data$Category == 'control_tool_hand_simple'] <- 4
  sorted_data$category_number[sorted_data$Category == 'control_tool_hand_complex'] <- 5
  sorted_data$category_number[sorted_data$Category == 'control_tool_face'] <- 6
  
  sorted_data$category_number[sorted_data$Category == 'meth'] <- 1 
  sorted_data$category_number[sorted_data$Category == 'meth_hand'] <- 2
  sorted_data$category_number[sorted_data$Category == 'meth_instrument'] <- 3
  sorted_data$category_number[sorted_data$Category == 'meth_instrument_hand'] <- 4
  sorted_data$category_number[sorted_data$Category == 'meth_injection_hand'] <- 5        
  sorted_data$category_number[sorted_data$Category == 'meth_face_activities'] <- 6
  
  #not using any meth_injection_instrument images
  sorted_data <- sorted_data[sorted_data$Category != 'meth_injection_instrument',]
  
  sorted_data$category_number[sorted_data$Category == 'opioid'] <- 1
  sorted_data$category_number[sorted_data$Category == 'opioid_hand'] <- 2
  sorted_data$category_number[sorted_data$Category == 'opioid_injection_instrument'] <- 3
  sorted_data$category_number[sorted_data$Category == 'opioid_instrument_hand'] <- 4
  sorted_data$category_number[sorted_data$Category == 'opioid_injection_hand'] <- 5
  sorted_data$category_number[sorted_data$Category == 'opioid_face_activities'] <- 6
  
  #number in each group that have been selected 
  number_selected <- c(0, 0, 0, 0, 0, 0)
  
  #go forward, marking the first 12 in each category of control images for selection
  #will set to TRUE for images selected to be used
  #this will be the 12 highest craving drug (opioid and meth separate) and lowest craving neutral images in each category
  sorted_data$selected <- FALSE
  for (i in 1:nrow(sorted_data)){
    #step forwards, taking the least craving-iducing neutral images
    if (sorted_data$ImageSet[i] != 'control') {
      next
    }
    this_category <- sorted_data$category_number[i]
    if (number_selected[this_category] > 11){
      next
    }
    number_selected[this_category] <- number_selected[this_category] + 1
    sorted_data[i, 'selected'] <- TRUE
  }
  
  
  
  #number in each group that have been selected 
  number_selected_meth <- c(0, 0, 0, 0, 0, 0)
  number_selected_opioid <- c(0, 0, 0, 0, 0, 0)
  #step through backwards, selecting most craving-inducing meth/opioid images
  #for (i in nrow(sorted_data):1){
  for (i in 1:nrow(sorted_data)){  
    
    this_category <- sorted_data$category_number[i]
    this_group <- sorted_data$ImageSet[i]
    if (this_group == 'meth'){
      if (number_selected_meth[this_category] > 11){
        next
      }
      number_selected_meth[this_category] <- number_selected_meth[this_category] + 1
      sorted_data[i, 'selected'] <- TRUE
    }
    if (this_group == 'opioid'){
      if (number_selected_opioid[this_category] > 11){
        next
      }
      number_selected_opioid[this_category] <- number_selected_opioid[this_category] + 1
      sorted_data[i, 'selected'] <- TRUE
    }
  }
  
  
  selected_data <- sorted_data[sorted_data$selected,]
  
  
  set.seed(12345678)
  selected_data$in_set <- NA
  #will randomize these
  in_set <- c(0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2)
  
  for (i in 1:6){
    selected_data[selected_data$category_number == i & selected_data$ImageSet == 'control', 'in_set'] <- in_set[sample(12)]
    selected_data[selected_data$category_number == i & selected_data$ImageSet == 'meth', 'in_set'] <- in_set[sample(12)]
    selected_data[selected_data$category_number == i & selected_data$ImageSet == 'opioid', 'in_set'] <- in_set[sample(12)]
  }
  
  library(tableone)
  
  allsets_tab <- CreateTableOne(vars = c('allsubjects_arousal_mean', 'allsubjects_craving_mean', 'allsubjects_typicality_mean', 'allsubjects_valence_mean', 
                                         'hue_mean', 'saturation_mean', 'value_mean'),
                                strata = c('ImageSet'), data = selected_data)
  print(allsets_tab)
  
  p1 <- as.numeric(substr(print(allsets_tab)[[30]], 2, 6))
  p2 <- as.numeric(substr(print(allsets_tab)[[31]], 2, 6))
  p3 <- as.numeric(substr(print(allsets_tab)[[32]], 2, 6))
  return(selected_data)
}

selected_data <- try_one_criteria(0.4, 0.6, 0.25,
                 0.3, 0.1, 0.5,
                 0.2, 0.2, 0.3)
##                                          Stratified by ImageSet
##                                           control      meth         opioid       p      test
##   n                                          72           72           72                   
##   allsubjects_arousal_mean (mean (sd))     2.53 (0.25)  4.72 (0.18)  4.65 (0.25) <0.001     
##   allsubjects_craving_mean (mean (sd))    11.88 (4.96) 88.54 (3.12) 88.39 (3.68) <0.001     
##   allsubjects_typicality_mean (mean (sd)) 12.22 (4.55) 89.70 (3.17) 89.54 (3.62) <0.001     
##   allsubjects_valence_mean (mean (sd))     4.89 (0.15)  4.42 (0.19)  4.44 (0.25) <0.001     
##   hue_mean (mean (sd))                     0.24 (0.18)  0.27 (0.14)  0.24 (0.13)  0.275     
##   saturation_mean (mean (sd))              0.26 (0.14)  0.30 (0.13)  0.27 (0.13)  0.139     
##   value_mean (mean (sd))                   0.42 (0.20)  0.39 (0.15)  0.41 (0.15)  0.574     
##                                          Stratified by ImageSet
##                                           control      meth         opioid       p      test
##   n                                          72           72           72                   
##   allsubjects_arousal_mean (mean (sd))     2.53 (0.25)  4.72 (0.18)  4.65 (0.25) <0.001     
##   allsubjects_craving_mean (mean (sd))    11.88 (4.96) 88.54 (3.12) 88.39 (3.68) <0.001     
##   allsubjects_typicality_mean (mean (sd)) 12.22 (4.55) 89.70 (3.17) 89.54 (3.62) <0.001     
##   allsubjects_valence_mean (mean (sd))     4.89 (0.15)  4.42 (0.19)  4.44 (0.25) <0.001     
##   hue_mean (mean (sd))                     0.24 (0.18)  0.27 (0.14)  0.24 (0.13)  0.275     
##   saturation_mean (mean (sd))              0.26 (0.14)  0.30 (0.13)  0.27 (0.13)  0.139     
##   value_mean (mean (sd))                   0.42 (0.20)  0.39 (0.15)  0.41 (0.15)  0.574     
##                                          Stratified by ImageSet
##                                           control      meth         opioid       p      test
##   n                                          72           72           72                   
##   allsubjects_arousal_mean (mean (sd))     2.53 (0.25)  4.72 (0.18)  4.65 (0.25) <0.001     
##   allsubjects_craving_mean (mean (sd))    11.88 (4.96) 88.54 (3.12) 88.39 (3.68) <0.001     
##   allsubjects_typicality_mean (mean (sd)) 12.22 (4.55) 89.70 (3.17) 89.54 (3.62) <0.001     
##   allsubjects_valence_mean (mean (sd))     4.89 (0.15)  4.42 (0.19)  4.44 (0.25) <0.001     
##   hue_mean (mean (sd))                     0.24 (0.18)  0.27 (0.14)  0.24 (0.13)  0.275     
##   saturation_mean (mean (sd))              0.26 (0.14)  0.30 (0.13)  0.27 (0.13)  0.139     
##   value_mean (mean (sd))                   0.42 (0.20)  0.39 (0.15)  0.41 (0.15)  0.574     
##                                          Stratified by ImageSet
##                                           control      meth         opioid       p      test
##   n                                          72           72           72                   
##   allsubjects_arousal_mean (mean (sd))     2.53 (0.25)  4.72 (0.18)  4.65 (0.25) <0.001     
##   allsubjects_craving_mean (mean (sd))    11.88 (4.96) 88.54 (3.12) 88.39 (3.68) <0.001     
##   allsubjects_typicality_mean (mean (sd)) 12.22 (4.55) 89.70 (3.17) 89.54 (3.62) <0.001     
##   allsubjects_valence_mean (mean (sd))     4.89 (0.15)  4.42 (0.19)  4.44 (0.25) <0.001     
##   hue_mean (mean (sd))                     0.24 (0.18)  0.27 (0.14)  0.24 (0.13)  0.275     
##   saturation_mean (mean (sd))              0.26 (0.14)  0.30 (0.13)  0.27 (0.13)  0.139     
##   value_mean (mean (sd))                   0.42 (0.20)  0.39 (0.15)  0.41 (0.15)  0.574
control_tab <- CreateTableOne(vars = c('allsubjects_arousal_mean', 'allsubjects_craving_mean', 'allsubjects_typicality_mean', 'allsubjects_valence_mean', 
                                       'hue_mean', 'saturation_mean', 'value_mean'),
                      strata = c('in_set'), data = selected_data[selected_data$ImageSet == 'control',])
control_tab
##                                          Stratified by in_set
##                                           0            1            2            p      test
##   n                                          24           24           24                   
##   allsubjects_arousal_mean (mean (sd))     2.50 (0.26)  2.56 (0.29)  2.54 (0.22)  0.739     
##   allsubjects_craving_mean (mean (sd))    11.17 (3.57) 11.61 (3.70) 12.84 (6.91)  0.488     
##   allsubjects_typicality_mean (mean (sd)) 12.06 (3.37) 11.69 (2.98) 12.92 (6.54)  0.636     
##   allsubjects_valence_mean (mean (sd))     4.87 (0.16)  4.90 (0.18)  4.90 (0.13)  0.775     
##   hue_mean (mean (sd))                     0.28 (0.19)  0.22 (0.18)  0.22 (0.14)  0.394     
##   saturation_mean (mean (sd))              0.26 (0.13)  0.25 (0.17)  0.26 (0.11)  0.966     
##   value_mean (mean (sd))                   0.42 (0.21)  0.39 (0.20)  0.45 (0.20)  0.595
meth_tab <- CreateTableOne(vars = c('allsubjects_arousal_mean', 'allsubjects_craving_mean', 'allsubjects_typicality_mean', 'allsubjects_valence_mean', 
                                    'hue_mean', 'saturation_mean', 'value_mean'),
                              strata = c('in_set'), data = selected_data[selected_data$ImageSet == 'meth',])
meth_tab
##                                          Stratified by in_set
##                                           0            1            2            p      test
##   n                                          24           24           24                   
##   allsubjects_arousal_mean (mean (sd))     4.71 (0.22)  4.74 (0.17)  4.70 (0.17)  0.715     
##   allsubjects_craving_mean (mean (sd))    87.35 (2.71) 89.26 (3.60) 89.00 (2.73)  0.070     
##   allsubjects_typicality_mean (mean (sd)) 88.34 (2.91) 90.30 (3.33) 90.46 (2.93)  0.034     
##   allsubjects_valence_mean (mean (sd))     4.42 (0.25)  4.40 (0.15)  4.44 (0.17)  0.817     
##   hue_mean (mean (sd))                     0.25 (0.12)  0.30 (0.18)  0.27 (0.12)  0.395     
##   saturation_mean (mean (sd))              0.29 (0.14)  0.30 (0.13)  0.32 (0.12)  0.683     
##   value_mean (mean (sd))                   0.36 (0.11)  0.43 (0.17)  0.38 (0.15)  0.204
opioid_tab <- CreateTableOne(vars = c('allsubjects_arousal_mean', 'allsubjects_craving_mean', 'allsubjects_typicality_mean', 'allsubjects_valence_mean', 
                                      'hue_mean', 'saturation_mean', 'value_mean'),
                           strata = c('in_set'), data = selected_data[selected_data$ImageSet == 'opioid',])
opioid_tab
##                                          Stratified by in_set
##                                           0            1            2            p      test
##   n                                          24           24           24                   
##   allsubjects_arousal_mean (mean (sd))     4.69 (0.23)  4.65 (0.25)  4.62 (0.26)  0.624     
##   allsubjects_craving_mean (mean (sd))    88.86 (3.29) 88.72 (4.17) 87.60 (3.56)  0.439     
##   allsubjects_typicality_mean (mean (sd)) 90.07 (2.99) 89.32 (4.52) 89.23 (3.26)  0.684     
##   allsubjects_valence_mean (mean (sd))     4.41 (0.23)  4.48 (0.27)  4.44 (0.27)  0.575     
##   hue_mean (mean (sd))                     0.23 (0.15)  0.24 (0.10)  0.24 (0.12)  0.865     
##   saturation_mean (mean (sd))              0.25 (0.12)  0.28 (0.13)  0.30 (0.15)  0.556     
##   value_mean (mean (sd))                   0.39 (0.13)  0.43 (0.17)  0.41 (0.16)  0.660
set0_tab <- CreateTableOne(vars = c('allsubjects_arousal_mean', 'allsubjects_craving_mean', 'allsubjects_typicality_mean', 'allsubjects_valence_mean', 
                                       'hue_mean', 'saturation_mean', 'value_mean'),
                      strata = c('ImageSet'), data = selected_data[selected_data$in_set == 0,])
print(set0_tab)
##                                          Stratified by ImageSet
##                                           control      meth         opioid       p      test
##   n                                          24           24           24                   
##   allsubjects_arousal_mean (mean (sd))     2.50 (0.26)  4.71 (0.22)  4.69 (0.23) <0.001     
##   allsubjects_craving_mean (mean (sd))    11.17 (3.57) 87.35 (2.71) 88.86 (3.29) <0.001     
##   allsubjects_typicality_mean (mean (sd)) 12.06 (3.37) 88.34 (2.91) 90.07 (2.99) <0.001     
##   allsubjects_valence_mean (mean (sd))     4.87 (0.16)  4.42 (0.25)  4.41 (0.23) <0.001     
##   hue_mean (mean (sd))                     0.28 (0.19)  0.25 (0.12)  0.23 (0.15)  0.515     
##   saturation_mean (mean (sd))              0.26 (0.13)  0.29 (0.14)  0.25 (0.12)  0.604     
##   value_mean (mean (sd))                   0.42 (0.21)  0.36 (0.11)  0.39 (0.13)  0.451
set1_tab <- CreateTableOne(vars = c('allsubjects_arousal_mean', 'allsubjects_craving_mean', 'allsubjects_typicality_mean', 'allsubjects_valence_mean', 
                                       'hue_mean', 'saturation_mean', 'value_mean'),
                      strata = c('ImageSet'), data = selected_data[selected_data$in_set == 1,])
print(set1_tab)
##                                          Stratified by ImageSet
##                                           control      meth         opioid       p      test
##   n                                          24           24           24                   
##   allsubjects_arousal_mean (mean (sd))     2.56 (0.29)  4.74 (0.17)  4.65 (0.25) <0.001     
##   allsubjects_craving_mean (mean (sd))    11.61 (3.70) 89.26 (3.60) 88.72 (4.17) <0.001     
##   allsubjects_typicality_mean (mean (sd)) 11.69 (2.98) 90.30 (3.33) 89.32 (4.52) <0.001     
##   allsubjects_valence_mean (mean (sd))     4.90 (0.18)  4.40 (0.15)  4.48 (0.27) <0.001     
##   hue_mean (mean (sd))                     0.22 (0.18)  0.30 (0.18)  0.24 (0.10)  0.186     
##   saturation_mean (mean (sd))              0.25 (0.17)  0.30 (0.13)  0.28 (0.13)  0.567     
##   value_mean (mean (sd))                   0.39 (0.20)  0.43 (0.17)  0.43 (0.17)  0.714
set2_tab <- CreateTableOne(vars = c('allsubjects_arousal_mean', 'allsubjects_craving_mean', 'allsubjects_typicality_mean', 'allsubjects_valence_mean', 
                                       'hue_mean', 'saturation_mean', 'value_mean'),
                      strata = c('ImageSet'), data = selected_data[selected_data$in_set == 2,])
print(set2_tab)
##                                          Stratified by ImageSet
##                                           control      meth         opioid       p      test
##   n                                          24           24           24                   
##   allsubjects_arousal_mean (mean (sd))     2.54 (0.22)  4.70 (0.17)  4.62 (0.26) <0.001     
##   allsubjects_craving_mean (mean (sd))    12.84 (6.91) 89.00 (2.73) 87.60 (3.56) <0.001     
##   allsubjects_typicality_mean (mean (sd)) 12.92 (6.54) 90.46 (2.93) 89.23 (3.26) <0.001     
##   allsubjects_valence_mean (mean (sd))     4.90 (0.13)  4.44 (0.17)  4.44 (0.27) <0.001     
##   hue_mean (mean (sd))                     0.22 (0.14)  0.27 (0.12)  0.24 (0.12)  0.402     
##   saturation_mean (mean (sd))              0.26 (0.11)  0.32 (0.12)  0.30 (0.15)  0.296     
##   value_mean (mean (sd))                   0.45 (0.20)  0.38 (0.15)  0.41 (0.16)  0.333
write.csv(selected_data, 'DCR_selected_split-MatchedHSV.csv', row.names=FALSE, quote = FALSE)

Craving vs RT scatterplot(s)

scatter_data <- subject_data[subject_data$RatingType == c('craving') & (subject_data$millisecondsOnPage < 20000) & subject_data$ImageSet %in% c('meth', 'opioid'), ]
ggplot(scatter_data, aes(x = Value, y = millisecondsOnPage / 1000, color = ImageSet)) + geom_point() + geom_smooth(method = 'lm') + 
  xlab('Craving') + ylab('Seconds') + 
  ggtitle(paste0('Craving vs Craving RT')) +
    scale_color_manual(values = cols)

scatter_data <- subject_data[subject_data$RatingType == c('craving') & (subject_data$ImageCompletionTime < 100) & subject_data$ImageSet %in% c('meth', 'opioid'), ]
ggplot(scatter_data, aes(x = Value, y = ImageCompletionTime, color = ImageSet)) + geom_point() + geom_smooth(method = 'lm') + 
  xlab('Craving') + ylab('Seconds') + 
  ggtitle(paste0('Craving vs Time On Image')) +
    scale_color_manual(values = cols)